Microscopic theory of the jamming transition of harmonic spheres 
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We develop a microscopic theory to analyze the phase behaviour and compute correlation func- 
tions of dense assemblies of soft repulsive particles both at finite temperature, as in colloidal mate- 
rials, and at vanishing temperature, a situation relevant for granular materials and emulsions. We 
use a mean-field statistical mechanical approach which combines elements of liquid state theory to 
replica calculations to obtain quantitative predictions for the location of phase boundaries, macro- 
scopic thermodynamic properties and microstructure of the system. We focus in particular on the 
derivation of scaling properties emerging in the vicinity of the jamming transition occurring at large 
density and zero temperature. The new predictions we obtain for pair correlation functions near 
contact are tested using computer simulations. Our work also clarifies the conceptual nature of the 
jamming transition, and its relation to the phenomenon of the glass transition observed in atomic 
liquids. 
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I. INTRODUCTION 



From atoms to balls 



It is common practice in physics textbooks to represent 
various states of matter by replacing atomic or molecu- 
lar arrangements by assemblies of spherical balls. The 
analogy goes much deeper, since simple theoretical mod- 
els of spherical, repulsive particles are known to display 
phase transitions between states of matter that are ac- 
tually observed in atomic systems. About 50 years ago, 
Bernal [1, 2] proposed a theory of the liquid state based 
on the intuition that dense liquids have a structure sim- 
ilar to dense assemblies of macroscopic grains, and per- 
formed with coworkers experiments to analyze the mi- 
crostructure of steel ball bearings. Since Bernal's insight- 
ful work, a large number of detailed experiments and nu- 
merical simulations have been made to characterize the 
structure of dense assemblies of spherical particles, both 
for hard particles relevant for granular materials [3] and 
for soft particles relevant for emulsions and foams [4]. 

However, while the structure of dense atomic liquids is 
by now well understood using the framework of statisti- 
cal mechanics [5], the problem of sphere packing proved 
much more difficult to solve [6, 7]. Indeed, there are 
two major obstacles to be faced when approaching this 
problem. First, sphere packings are inherently nonequi- 
librium systems where thermal fluctuations play no role. 
This has for instance led researchers to extend concepts 
of statistical mechanics to specifically address athermal 
packings implementing ideas first put forward in partic- 
ular by Edwards [8-11]. A second major difficulty is that 
one must deal with both fluid and solid disordered states 
of matter, which are poorly understood subjects even for 
atomic systems driven by thermal fluctuations. In partic- 



ular, while dense liquids are appropriately described by 
liquid state theory, viscous liquids and amorphous mate- 
rials are still the focus of intense research [12]. Following 
Bernal's insight, it was suggested that the fluid to solid 
transition observed in athermal systems shares similar- 
ities with the glass transition of supercooled molecular 
fluids [13]. Since its publication, the "jamming phase di- 
agram" has striked lively debates and a large literature 
seeking in particular to understand similarities and dif- 
ferences in parallel studies of both thermal and athermal 
amorphous media [14]. 

In this work, we are interested in describing theoreti- 
cally fluid and solid disordered states of matter observed 
in dense assemblies of spherical repulsive particles both 
for finite temperatures relevant for colloidal materials, 
and in the athermal limit relevant for emulsions and 
foams. Unlike Bernal, we now have at our disposal theo- 
ries for the liquid state, but this is not sufficient [15, 16]. 
Liquid state theory is essentially based on perturbation 
theory extrapolating from the gas state, and for this 
reason it inherently fails to describe the phase transi- 
tions to low temperature and high density solid phases 
of matter. However, the liquid-crystal transition is con- 
veniently described by density functional theory, where 
one starts from an approximated form for the free en- 
ergy as functional of the local density profile and looks 
for a transition between uniform and periodically mod- 
ulated minima of this functional. Early studies of the 
glass transition [17, 18] used density functional theory 
and looked for amorphous density profiles that minimize 
the free energy. Theoretical developments in the field of 
spin glasses [19-21] suggested similarities between super- 
cooled liquids and a class of mean-field spin glass models 
which can be extensively analyzed theoretically [22, 23]. 
These works gave birth to the so-called Random First 
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Order Transition (RFOT) theory [21], which is currently 
one of the most consistent theories of the glass transi- 
tion [12, 24-26]. Along with these conceptual develop- 
ments, RFOT theory has been turned into a quantitative 
computational scheme in a series of pioneering works by 
Monasson [27], Franz, Mezard and Parisi [28-31]. These 
works can schematically be seen as an extension of the 
tools developed in liquid state theory to study systems 
near the glass transition and the properties of amor- 
phous solids, but they remain "mean- field" in nature in 
the sense that the glassy phase is assumed to be of the 
same nature as the one obtained from various mean-field 
perspectives [12]. We wish to emphasize, however, that 
"mean-field" has a very specific meaning in the context 
of glass theories. Although the theory we shall develop 
is a mean-field theory, it remains fully microscopic in na- 
ture (wc start from the interaction between the particles), 
and it is able to capture many-body correlations between 
particles giving rise to a complex free energy landscape, 
which is crucial to accurately predict both the phase dia- 
gram of the system and correlation functions. In particu- 
lar, we are able to provide a theoretical description of the 
distinct nature of both glass and jamming transitions. 

Our goal is thus to apply concepts and tools from 
RFOT theory to analyze the glassy and jammed phases 
of systems of soft repulsive particles, thus contributing 
to better understanding amorphous, athermal states of 
matter [14, 32]. In previous work, Parisi and Zamponi 
(PZ) extended the scheme devised by Mezard and Parisi 
(MP) to treat systems of hard spheres [33]. The result is a 
RFOT theory of the glass transition and jamming of hard 
spheres; this theory is expected to be correct in the limit 
of large space dimension [33] and for a modified model 
where particles are randomly displaced over the box [34] 
(both limits being of a mean- field nature), and gives 
quantitatively correct results in three dimensions [33]. 
In a similar spirit, Mari et al. devised a model of hard 
spheres with diluted interactions [35] where "mean-field" 
replica calculations become exact, although the solution 
remains quite involved [36]. We wish to unify the ap- 
proaches of MP and PZ to describe simultaneously the 
critical properties on both sides of the jamming transi- 
tion occurring in athermal materials, and in its vicinity 
at finite temperatures. The theory that we develop in 
this work is, in our opinion, the equivalent, for the jam- 
ming transition, of a Landau mean-field theory for second 
order phase transition. 



B. Glass and jamming transitions in harmonic 
spheres 



A model system that can interpolate between finite 
temperature glasses and hard spheres is the model of har- 
monic spheres. An assembly of N spherical particles of 
diameter a is enclosed in a volume V in three spatial di- 
mensions, interacting with a soft harmonic repulsion of 




FIG. 1: The jamming transition of harmonic splieres at zero 
temperature. Below ipj, particles do not overlap (as in hard 
sphere configurations) and the system flows. Above tpj , there 
is a finite density of overlaps, and the energy and pressure are 
finite. In this region, hard spheres configurations cannot be 
found. 



finite range: 



(r) = e(l-r/cr)2 6l((T-r) 



(1) 



where r > is the interparticle distance, 9{r) is the 
Heaviside step function and e controls the strength of 
the repulsion. The model (1), originally proposed to de- 
scribe wet foams [37], has become a paradigm in numer- 
ical studies of the T = jamming transition [32, 38]. It 
has also been studied at finite temperatures [39-41], and 
finds experimental realizations in emulsions, soft colloids 
and grains [42]. The model has the two needed control 
parameters to explore the jamming phase diagram: the 
temperature T (expressed in units of e), and the fraction 
of the volume occupied by the particles in the absence of 
overlap, the volume fraction Lp = ixNa'^ / {QV). Wc now 
set a and e to unity. 

Over the last decade, a large number of numerical ob- 
servations was reported for this model [14, 32]. A jam- 
ming transition is observed at T = at some critical 
volume fraction (pj^ the density above which the pack- 
ings carry a finite density of particle overlaps. This tran- 
sition is pictorially represented in Fig. 1. Numerically, 
the zero-temperature energy density, cqSj and pressure, 
P, are found to increase continuously from zero above 
(pj as power laws [38]. The pair correlation function of 
the density fluctuations [5], g{r), develops singularities 
near ipj [43, 44], which are smoothed by thermal exci- 
tations [41]. In particular, 17(1) = 00 at ipj and T = Q. 
This behavior implies that the density of contacts be- 
tween particles, z, jumps discontinuously from to a fi- 
nite value, Zj, at ipij. and increases further algebraically 
above (pj [32, 38]. Thus, the jamming transition ap- 
pears as a phase transition taking place in the absence 
of thermal motion, with a very peculiar critical behav- 
ior and physical consequences observable experimentally 
in a large number of materials. Note however that in 
the unjammed state below (pj, the specific configurations 
sampled by the system might depend on the details of 
the preparation protocol. In the absence of fluctuations 
for instance [45, 46], the contact number is not vanishing 
even below the jamming point, because particles remain 
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FIG. 2: (Color online) A schematic phase diagram for the 
glassy states of the model defined by Eq. (1) showing both the 
glass transition line at Txif) and the jamming point at T = 
and ifi — ifiGCP- The Parisi-Zamponi (PZ) approach only 
holds at T = and ifi < ipacp , while the Mezard-Parisi (MP) 
approach holds for large enough densities and temperatures. 
In this paper we develop a method to fill the gap between the 
PZ and MP analytical schemes to explore the vicinity of the 
jamming transition. 



in contact and the left cartoon in Fig. 1 is not a faithful 
representation. 

As discussed in detail by Krzakala and Kurchan [47], 
the jamming transition is very similar in nature to the 
so-called SAT-UNSAT transition of random constraint 
satisfaction problems [48, 49]. These are problems where 
a large set of constraints has to be simultaneously sat- 
isfied. Upon increasing the density of constraints per 
variable, a phase transition is found, from a phase where 
the problem is typically solvable (SAT) to a phase where 
it is typically not (UNSAT). In the present context, the 
constraints are the non-overlapping conditions between 
spheres, and the control parameter is the density. At ipj, 
the system goes from a phase where the non-overlap con- 
ditions can be satisfied to a phase where they cannot. In 
the context of random constraint satisfaction problems, 
it is very useful to introduce a soft version of the con- 
straints, study the problem at finite temperature T, and 
finally take the T — >■ limit [50-52]. The reason is that 
introducing a temperature allows to use powerful statis- 
tical mechanics tools, that can be used in a context where 
they are not a priori relevant [47]. 

To address the purely geometrical packing problem of 
soft spheres, we then suggest to study first the statisti- 
cal mechanics of the model defined by Eq. (1) at finite 
temperatures, before taking the T — >■ limit where jam- 
ming occurs. Based on the analogy with constraint sat- 
isfaction problems, we have in mind the schematic phase 
diagram reported in Fig. 2: see in particular the result 
for Potts glasses in [52], and see also [53]. In this pic- 
ture, that we will derive in the following, the liquid un- 
dergoes a glass transition at some temperature Tk{^)- 
The line Tk{^) vanishes at a volume fraction which 
corresponds also to the ideal glass transition for hard 



spheres [33] . The point Lpx is not the jamming transition: 
indeed, the ground state energy and pressure remain zero 
across fK- Above i^k, the system enters at zero temper- 
ature a hard sphere glassy state. In this state, particles 
vibrate near well-defined (but random) positions, and the 
system is not yet jammed [33]. Jamming happens when 
the glass reaches its close packing density, which we shall 
call "glass close packing" (GCP) [33]. We identify the 
jamming transition with (pccp- In Fig- 2 the density 
interval [^Kt'Pgcp] is simply the amorphous analog of 
the interval (p G [0.54,0.74] for ordered states of hard 
spheres where a compressible crystalline structure exists 
at thermal equilibrium. 

The main theoretical difficulty in the study of the jam- 
ming transition at (pacp is that it happens deep inside 
the glass phase. Therefore, we need to develop first an 
accurate theoretical description of the glass phase. As 
discussed in detail in Ref. [16], previous theories of the 
glass phase fail (see Fig. 2). The MP approach holds only 
at high enough temperature or density, where a simple 
harmonic approximation for vibrations in the glass holds. 
The PZ approach holds only in the hard sphere region, 
T ~ and ip < (fccp- The central theoretical achieve- 
ment that we report in this paper is an approximation 
scheme that can naturally interpolate between those of 
MP and PZ and is correct in the entire vicinity of the 
jamming transition, thus allowing us to fully explore the 
phase diagram of Fig. 2, and in particular the region of 
low temperature, T ^ Tk{ip) and tp ~ focp- We ob- 
tain theoretically from a first principle calculation a large 
number of the observed behaviors of harmonic spheres, 
and predict new results for the correlation functions of 
this system around the jamming transition. A short re- 
port of our results has been published in Ref. [54] . 



C. Discussion 

Our approach is very different from, but complemen- 
tary to, recent theoretical works on jamming. 

The aim of our work is to show, directly from the 
Hamiltonian, that the jamming transition exists, and 
determine from first principles its location and proper- 
ties. Other approaches assume the existence of jammed 
states and try to obtain geometrical informations on 
them [10, 55]. A particularly interesting approach is 
to develop a scaling picture of the jamming transition 
by showing that the transition is characterized by a di- 
verging correlation length [56-58] . This approach allows 
to obtain a very detailed description of the jamming 
transition [59], that explains the anomalous scalings in 
shear [38, 60] and transport [61] properties of amorphous 
jammed packings. 

Moreover, our method is crucially based on equilib- 
rium statistical mechanics. In the glassy phase, on ap- 
proaching Tk{p) and below it, the equilibrium Gibbs 
distribution of the system is the superposition of many 
distinct ergodic components, each representing a differ- 
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ent glassy state; each of these "basins" (see Appendix A 
for a more detailed discussion) is visited, in our theory, 
according to its equilibrium weight. On the contrary, 
most of the numerical works on jamming focus on specific 
dynamic protocols, that often are athermal and there- 
fore always out of equilibrium. It has been shown that 
these protocols sample jammed states with highly non- 
uniform probabilities [46, 62]. Even thermal protocols 
(see e.g. [43]) inevitably fall out of equilibrium before 
Tx{^) is reached, and therefore produce jammed pack- 
ings with probabilities that are strongly different from 
the equilibrium ones [33, 47]. Still, a crucial observation 
is that the critical properties around the (protocol de- 
pendent) jamming density ipj are largely independent of 
the protocol, and therefore of the probability with which 
jammed states are visited. This observation suggests that 
the critical properties are shared by almost all jammed 
states, and therefore an equilibrium sampling of these 
states can appropriately describe these properties, as we 
showed in this work. Of course, some care will have to 
be taken when comparing with some specific protocols, in 
which no fluctuations are present at all (neither thermal, 
nor mechanical) [45, 46]. In those cases, as we already 
mentioned, the contact number is not vanishing even be- 
low the jamming point. Once these effects are taken into 
account, we believe that the generic picture we discuss in 
the following also applies to these protocols. 

Finally, it is worth mentioning that if compression 
rates are small, partial crystallization takes place and the 
system can jam at any packing fraction between ipccp 
and the close packing density corresponding to the most 
dense crystal [63-65]. In our theory, crystallization is 
not included by construction, as we will explain in the 
following, therefore this effect is absent. Extending this 
study to take into account partial crystallization would 
be interesting but probably very hard. 

The outline of the paper is the following. In Sec. II, we 
present the approximation scheme developed earlier, and 
describe our own approximation specifically developed to 
study harmonic spheres. In Sees. Ill, IV and V respec- 
tively, we present results concerning the glass transition, 
the jamming transition and the correlation functions. In 
Sec. VI we compare the theoretical results with numerical 
data. We discuss and conclude the article in Sec. VII. 



sitions upon lowering the temperature towards the glass 
state [23]: at a first transition Td, the metastable states 
appear, giving rise to a finite complexity S = N^^ log A/", 
which is the entropic contribution due to the exponen- 
tial number of metastable states. Below this tempera- 
ture, the dynamics of the system is arrested: this tem- 
perature is indeed identified with the dynamical transi- 
tion predicted by Mode-Coupling Theory [20, 25]. At a 
lower temperature Tk, the complexity vanishes, and the 
number of metastable states is sub-exponential; this has 
been identified with the Kauzmann entropy crisis transi- 
tion [66]. Therefore at the mean-field level the complex- 
ity is finite between Td and Tk- Although we will not 
discuss this physics in the rest of the paper, it is worth 
noting at this point that RFOT theory predicts that in 
a three (or finite) dimensional system the complexity is 
finite only locally. This is because the system gains en- 
tropy by accessing many different states; this creates an 
entropic driving force that allows the system to main- 
tain its ergodicity by visiting all its metastable states. 
A nucleation theory of supercooled liquids has been de- 
veloped based on these ideas [24-26] , which predicts that 
the transition at Td is avoided in finite dimension through 
activated events, and the dynamics is arrested only at 
Tk. These ideas have been initially proposed on a phe- 
nomenological basis [24-26] but are now starting to be 
confirmed by renormalization group computations in sim- 
ple models [67, 68] . While nucleation processes are surely 
important to describe activated dynamics below Td, one 
expects that the local structural properties will not be 
affected by them and therefore mean-field theory should 
provide precise predictions as long as static and local ob- 
servables are concerned. For this reason, in this work we 
follow Mezard and Parisi [28, 30] and directly apply the 
mean-field concepts described above to three dimensional 
glass-formers, neglecting activation. 

Let us denote by C the microscopic configuration of 
the system. The starting point of the discussion is the 
assumption that in the glassy region, phase space can be 
partitioned in pure states, indexed by a = 1 , • • • , A/", in 
such a way that a given configuration C belongs to one 
state a. A practical implementation of this construction, 
which has been studied recently in Ref. [69] , is discussed 
in Appendix A. Under this decomposition, the partition 
function can be written as: 



II. THE METHOD 
A. The complexity 



^^^g-/3£;(c)^^^g-mc)^ (2) 

C Q=lC6a 



In mean-field models of glasses, upon decreasing tem- 
perature, a number Af of metastable states, that scales 
as the exponential of TV, appears [23]. The dynamics of 
the glass is slowed down considerably: once the glass has 
fallen into one of these metastable states, the barrier it 
has to cross in order to jump to another metastable state 
is proportional to N, due to the mean-field nature of the 
models. In mean-field theory there are two distinct tran- 



We define an effective free-energy of a pure state /„ 
by fa = -i7;aln(I]ceae~'^^^'^^)' and the complexity 
S as the logarithm of the number J\f{f) of pure states 
at a given free-energy /: S(/) = ■^logA/'(/). We 
assume that, as in exactly solvable mean-field mod- 
els [23], the complexity S(/) is non-zero in a finite in- 
terval [frnin, fmax] of frcc energy, and it vanishes linearly 
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at frain- Wc find: 

Z= f dfN{f)e-Pf = I d/e-^^(/-^^(^)). (3) 



It is convenient for the following of the paper to consider, 
instead of the total free energy F ~ ~TN~^ log Z , a 
"free entropy" S = ^ogZ = —f3F; the advantage is 
that this quantity has a finite limit for T — ^ 0, which 
is the entropy of hard spheres. Saddle-point evaluation 
at thermodynamic limit gives an expression for the free 
entropy: 

S{T, ^) ^ 1 In Z = S](/*(r, ^)) - /3/*(T, (4) 



|^(r(T,^)) = i. 



(5) 



The glass transition is reached at the Kauzmann temper- 
ature Tk{lp), when the complexity vanishes: 



= 0. 



(6) 



When T < TK{f), the system has the free energy 
fmin(T^ if) and the complexity sticks to zero. There- 
fore, detecting the glass transition and calculating the 
free energy of the glass phase amounts to finding a way 
to compute the complexity as a function of the free en- 
ergy. 

A schematic plot of Tk{^p) for harmonic spheres is re- 
ported in Fig. 2. Although our approach does not allow 
for a reliable calculation of Td{(p) [33], we expect that 
the latter will be slightly bigger than Tk{^) and follow 
a similar trend, as it can be shown in mean-field lattice 
models (see e.g. Fig. 2 in Ref. [70]). 



B. Replicating the system 

The need to replicate the system in order to compute 
the thermodynamics of the glass originates from the dif- 
ficulty in defining an order parameter for the glass tran- 
sition. Indeed, below Tk the system freezes in an amor- 
phous density profile which is structurally very similar to 
a typical liquid configuration, therefore distinguishing the 
liquid and the glass by means of some observable related 
to structure (like a correlation function) is extremely dif- 
ficult. Edwards and Anderson proposed a way out of 
this difficulty by introducing an order parameter defined 
as [22] 



(7) 



where {6p{x)) — {p{x)) — p is the thermodynamic average 
of the fluctuation of the local density p{x) with respect 
to its spatial average p = N/V. This order parameter 
is clearly zero in the liquid phase while it is non-zero 
in a frozen amorphous phase. The problem is that its 
computation is not easy since one has to compute the 



local density and average its square. A solution is to 
introduce two identical copies of the system and compute 



1 
V 



d?x{8p^{x)){5p2{x)) 



(8) 



The problem is that if the two replicas are independent, 
they will freeze in completely different profiles, and qsA 
will be zero since the integrand will have random signs. 
To avoid this problem, we must couple the two replicas 
via a small interaction, and then compute qea by sending 
first iV — > oo and then switching off the interaction. In a 
liquid phase, since ergodicity is maintained, the copy will 
eventually decorrelate from the original, and the order 
parameter will go to zero. In an ergodicity-broken phase, 
the copy will stay trapped near the original and the order 
parameter will go to a finite value. This prescription 
allows one to detect the glass transition but does not 
give any information on the free-energy of the glass and 
thus the glass phase itself. 

Monasson [27] showed that replicas can provide a sim- 
ple way to compute the complexity. He introduced an ar- 
bitrary number of replicas to, and assumed that a small 
coupling between each pair of replicas is present, in such 
a way that all replicas are constrained to be in the same 
metastable state. Since the number of metastable states 
does not change, the complexity is unaffected by this 
procedure. The partition function Z,„ of this replicated 
liquid then reads: 



dfe 



-NP{mf-Ti:(f)) 



(9) 



The number to of replicas is eventually taken to be non- 
integer, creating an additional to dependence of the free- 
energy and complexity of the glass. Equations (4) and 
(5) become now: 



5(to; T, ^) = nr ("^; T, ^)) - mprim; T, (10) 

(11) 



|^(r(TO;T,^)) = ^. 



This TO-dependence of the free entropy allows one to com- 
pute the complexity as follows. We can determine E and 
/* as functions of to at fixed (ip, T) using 



S = -TO- 



:d{S/m) 



dm 



dm 



(12) 



Then one can eliminate parametrically the dependence 
on TO and in this way reconstruct the full curve S(/). 

Wc showed in the previous section that for T < Tk 
the free energy of the glass is equal to fmin {T) , the point 
where S(/) = 0. Therefore, to compute fminiT), we 
must impose that S(/) — ^^g^-* = 0, which is equivalent 
to finding the minimum of 5/to with respect to to. Let us 
call TO* (T, if) the point where S/m assumes its minimum. 
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Since E = at the minimum, using Eq. (10), we get 
S{m*;T,ip) = —'m*l3fmin{T,'p) which imphes that: 



Sglass{T, (fi) 



-l3fglassiT,(p) = - 

S{m*;T,(p) 

— r— = miE 

m*[T,ip) m 



-i3fmin{T, if) 

S{m;T, (p) 



(13) 



Now we see that, in order to compute the thermodynam- 
ics of the ideal glass, one must be able to calculate the 
free entropy of a replicated liquid. In the following, we 
will present several approximation schemes that we ap- 
plied and adapted, when necessary, for harmonic spheres. 



C. Effective potential approximation 

Computing the free-energy and complexity of the 
glass amounts to computing the free entropy of a m- 
times replicated liquid. One can use a replicated ver- 
sion of standard liquid theory approximations such as 
the replicated Hyper-Netted-Chain (HNC) approxima- 
tion [28, 29], but these kind of approximations usually 
break down deep in the glass phase. For a discussion of 
the reasons for the failure of replicated HNC, see Ref. [33]. 
One can also follow Refs. [30, 31] and start at high den- 
sity, and perform cage expansions: at high densities the 
copies stay close to the originals, forming molecules of 
size A. One can then expand the replicated free entropy 
with respect to this parameter A. However, to be able 
to explicitly compute the various integrals appearing in 
the computation, one traditionally has to suppose that 
the attraction between copies is harmonic, an assump- 
tion which breaks down for hard spheres, as discussed 
in [33]. To treat correctly the replicated hard sphere 
system, Parisi and Zamponi [33] developed an effective 
potential approximation, which in the end amounts to 
performing an expansion in powers of y/A instead of A. 
We explain in the following the extension of the effective 
potential approximation to finite temperature harmonic 
spheres. By construction, we will see that we are bound 
to recover both high density MP results and zero tem- 
perature PZ results in a unified treatment. 

The starting point of all these approaches is the as- 
sumption that (due to the implicit coupling between 
replicas) the replicated system is composed of molecules 
made of m atoms, each belonging to a different replica. 
Making use of standard theory of molecular liquids, one 
can express the free energy as a functional of the single 
molecule density and the interaction potential, see [33] 
for details. In order to make the computation tractable, 
one then makes a Gaussian ansatz for the probability dis- 
tribution function of the positions x of the replicas within 
a molecule [31, 33]: 



pixi ■■■x^)= / (fx Y[ 



e 2X- 



(14) 



We can now integrate out all replicas except one. The ef- 
fective interaction between particles in replica 1, 4'eff, is 




FIG. 3: (Color online) A schematic representation of the ef- 
fective potential approximation. Each particle in the original 
liquid is replicated m times (dashed spheres). Assuming that 
the replicated particles form a molecule of average cage size 
A, we trace out in the partition sum the degrees of freedom 
of (m — 1) copies of the liquid to obtain an effective one- 
component liquid (black spheres) interacting with an effective 
pair potential 0e// (green lines). 



obtained by averaging the full interaction over the prob- 
ability distribution of the two molecules: 



\a=2 



(fx2 ■ ■ ■ d^Xmcfy2 ■ ■ ■ d^Vm'X 
ni 

X p(fi • • • x„,)p{m ■ ■ ■ y™) n e-'3^(^--y~") 



o=l 



(15) 



We get after some simple manipulation: 



(16) 

where j2A is a normalized and centered Caussian of vari- 
ance 2A and g(A,T;r) = / d^r'-f2A{r^)e-^'>''^^~^^ is a 
function that can be explicitly computed for our har- 
monic potential and is given in Eq. (B2) of Appendix B 1. 
This calculation is schematically represented in Fig. 3. 
We can rewrite this result using bipolar coordinates: 



X / du u 

10 



:X 



ryAnA 



q{A,T;u) 



m— 1 



(17) 



The integration of replicas 2 • • • m induces also three- 
body (and more generally many-body) interactions on 
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replica 1; if we neglect these, and just keep the two-body 
effective potential, we obtain that the free entropy of the 
replicated liquid is the sum of an harmonic part 5/i(m, A) 
plus the free energy of a liquid interacting via the effective 
potential: S = Sh + Suqlipeff] [33]. Furthermore, we 
define 



Q{r) — e^'^^'^o-fJ'*'')"'"'^*'')' — 1 



(18) 



we note that Q{r) is equal to zero when A = 0, therefore 
for small A we suppose that Q{r) is a small perturba- 
tion, and wc make use of equilibrium liquid perturbation 
theory. Wc finally obtain for the replicated free entropy 
S{m,A;T,(p) = ^logZ„: 

S{m,A;T,ip) Sh{m,A) +Suq{T/m,ip) 



d-'rguq{T/m,ip;r)Q{r) , (19) 



3 3 
Sh{m,A) = -(m - l)ln(27rA) + -(m - 1 + Inm) 



(20) 



where Siiq{T, cp) = N^^ log Zuq is the free entropy of the 
liquid, giiq(T,ip;r) its correlation function (to be com- 
puted at temperature T/m). 

Equation (19) allows for a full calculation of the repli- 
cated free entropy, and gives back both the MP small cage 
expansion in powers of A at high density, and the PZ ex- 
pansion in powers of y/A at zero temperature. Details 
about this can be found in the Appendix B 2 and C. The 
main theoretical goal of the paper is now achieved [54] : 
we obtained a set of equations that allow us to make a 
connection between zero temperature hard spheres and fi- 
nite temperature harmonic spheres. Note that the molec- 
ular liquid is assumed to be in a liquid state, therefore 
crystallization is excluded a priori from the theory [33]. 



D. Low temperature approximation 

Our purpose in interpolating between zero and finite 
temperature was to be able to probe the vicinity of the 
jamming point in the (i^, T) phase diagram, thus we are 
interested mainly in the low temperature behaviour of 
the glass. We can exploit this by making an approxima- 
tion that will allow us to push the analytical calculation 
much further. Defining the cavity distribution function 
y by g{T, (p; r) = e~^'^^'^'^y{T, (p; r), wc make the approx- 
imation of taking the cavity function of the liquid as a 
constant and equal to its T = 0,r = 1 value, which we 
cah yl^qif)- As T goes to zero, we see that the expo- 
nential factor converges towards a step function around 
r = 1, so that in all integrals that are cut at r = 1 by the 
pair potential, we can safely evaluate y at its r = 1 value. 
Expanding y in powers of T, it is easy to convince oneself 
that the temperature dependence leads to subdominant 



contributions in the integrals. Thus we suppose: 

guq{T/m, v.; r) ^ e"^'"*W2/,f/(^). (21) 

Plugging this approximation into Eq. (19), we get (the 
details of the computation are in Appendix B 3): 

S{m,A;T,(p) = Shim, A) + Siiq{T/m,ip) 

+ 4^ySq^{^)G{m,A;T), (22) 

G{m,A;T)^3 drr^[q{A,T;r)"' - e-^''"^^% (23) 



where Suq{T/m,Lp) is the free entropy of the non- 
replicated liquid, evaluated at an effective tempera- 
ture T/m. These are the finite-temperature version of 
the replicated free entropy of hard-spheres obtained in 
Ref. [33]. 

Many technical aspects presented in the following al- 
ready appear in [33], apart from the existence of a new 
control parameter, the temperature. The average cage 
radius A*(m;T, tp) is obtained by maximization of S 
against A. We find an implicit equation for A*{rn]T,Lp): 



J{m, A* (m; T, ip)] T) = 
J{ra,A]T) 
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A dG{A,m;T) 



1 — TO dA 



1 — m 



drr^q{A,T;ry 



(24) 

(25) 

,_, dq{A,T;r) 
dA 



E. Liquid theory 

The last quantities that we need to compute are the 
free entropy Suq and the pair correlation function g^q (r) 
of the liquid. Given the pair correlation function, which 
represents the probability of finding a particle at distance 
r from a particle fixed at the origin, we can express the 
internal energy U as: 

/>oo 

U{T,(p)^l2(p drr^g{T,(p;r)(P{r) . (26) 
Jo 

Plugging the low temperature approximation Eq. (21) for 
the liquid into this, we obtain: 

Uuq{T,^)^^12^y"'{^) C drr\l-rfe-P<^'-^^' . 



(27) 



We are interested in computing the free entropy Suq of 
the liquid, which is related to the liquid free-energy Fuq 
by Siiq{T, ip) = —f3Fiiq(T, ip). The free entropy has a fi- 
nite limit for hard spheres T — >■ 0, thus allowing to make 
the connection between hard spheres and soft spheres at 
finite temperature. Making use of the standard identity 



8 



U{T) = ^ixTTJ' '^^ derive the low temperature ap- 
proximation for the free entropy of the hquid: 



^VT(2 + T)crf^^^^ 



-T 



-l/T 



, (28) 



where yflqi'-p) and S[^^{(p) are short-hand notations for 
yiiq{T = 0, v?; r = 1) and SuqiT = 0,ip). 

At this level of approximation, it is clear that the only 
input that is needed from liquid theory is the equation of 
state of the hard sphere liquid. From any given equation 
of state one can easily deduce the hard sphere free en- 
tropy 5;^"^ and cavity function yf^^ ■ The most reasonable 
choice would be to use the phenomenological Carnahan- 
Starling (CS) equation of state, as in [33], that provides 
the best fit to the hard sphere pressure. 

In this paper we use instead the Hyper-Netted Chain 
(HNC) approximation: the HNC equations are non-linear 
integro-difFcrcntial equations that require a numerical so- 
lution. Although HNC is known to be less accurate for 
the hard sphere system, using HNC allows us to also com- 
pare the present directly our results to the ones obtained 
from the MP approach valid at large density and finite 
temperatures. The procedure for solving these equations 
is described in full details in [71]. The HNC approxima- 
tion overestimates yff^ by 20% in the relevant range of 
volume fraction (p ~ 0.64. This has the effect of reduc- 
ing the glass transition density obtained from the the- 
ory, from the value (pK = 0.62 obtained from CS [33] 
to (pK = 0.58 obtained with HNC. Therefore the reader 
should keep in mind that the glass densities reported in 
the following are lower than the correct ones. In any case, 
here we are more interested in the low-temperature scal- 
ing in the glass phase than to the actual value of the glass 
transition density. We also checked that the scaling re- 
sults are insensitive to the precise choice of the equation 
of state of the hard sphere liquid. 



III. THERMODYNAMICS OF THE GLASS 



A. Complexity and phase diagram 



As discussed above, the glass transition is signalled by 
the point where the saddle point in Eq. (3) reaches the 
minimum fmin at which the complexity vanishes. In 
the replica formalism, the equilibrium complexity 
in Eq. (6) corresponds to the complexity in Eq. (12) eval- 
uated at m = 1. We call this quantity the "equilibrium" 
complexity of the liquid Eeq(T, (p). The latter is easily 
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2x10-5 4xl0->"^ 6x10-5 8x10-5 ixlO"^ 
T 



FIG. 4: (Color online) Equilibrium complexity against tem- 
perature for several volume fractions around ifK, calculated 
within the elfective potential approximation. Inset: complex- 
ity of hard spheres, obtained as the T — >■ limit of the com- 
plexity Eeq, plotted against the volume fraction 



computed by expanding the equations around m ~ 1: 



I]e,(T, ^) = Sug{T, (p) - - ln(2^A*(m = 1: T, ^) - 3 
- 12ipy"^{ip)H{m = l,A*{m = l;T,(p);T), 



(29) 



H{m,A-T) 



1 dG{m,A-T) 
m dm 

/•OO 



= — / drr2g(A,r;r)'"lng(A,r;r) 
"Wo 

_ A /V2^(r)e-"''*W. (30) 
Jo 

We see that the only unknown in the previous equation is 
the optimal cage radius A* . It is computed with Eq. (24) 
evaluated at m = 1. 

Since our computation makes connection with the hard 
spheres calculation of PZ, we expect to retrieve all their 
results in the T = limit for ip < ipccp- Thus for 
ip smaller than pK (the Kauzmann transition of hard 
spheres), the system should be in a liquid phase, and thus 
its complexity should not vanish, and converge to the 
hard sphere complexity Yi^^^ (p) found by PZ. In practice, 

since PZ performed a -s/A expansion, whereas we did not, 
our numerical values for (pK and T,^^ {p) could differ 
slightly. Nevertheless, A is typically found between 10"'^ 
and 10~^, so that the difference is indeed very small. For 
example px calculated by PZ differs from our value by 
less than 6 • 10"^. 

For p > ipK , the complexity vanishes at a temperature 
Tk{'-p), called the Kauzmann temperature, that increases 
with the density. These results are summarized in Fig. 4, 
where we show the complexity as a function of tempera- 
ture for several volume fractions around (pK ■ In the inset 
we show the zero temperature limit of the complexity 
S^'^ as a function of the density, that vanishes at p}K- 
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FIG. 5: (Color online) Kauzmann temperature against vol- 
ume fraction, within the the effective potential approach de- 
veloped in this work (open squares). The quadratic fit (red 
hne) gives an estimated Kauzmann transition for hard spheres 
at i^K ~ 0.576898. The small cage Mezard-Parisi approach 
breaks down at low temperatures and low densities, as ex- 
plained in Ref. [16]. 



From the complexity we deduce the phase diagram 
shown in Fig. 5. In this figure we report Tfj-((p), as ob- 
tained in the framework of the effective potential approx- 
imation developed in this paper, and we compare it with 
the result obtained from the Mczard-Parisi small cage ex- 
pansion [16]. In the effective potential case we find that 
the Kauzmann temperature goes to zero at ipx ~ 0.5769, 
quadratically in 1^9 — ipx- Conversely, the result from the 
small cage expansion is a finite Tn^ip) which jumps to 
zero abruptly at a value of if which is unrelated to hard 
sphere results. This is due to fact that the small cage 
expansion is valid only in the region indicated in Fig. 2. 
On the other hand, our effective potential computation 
becomes inaccurate when the temperature is too high be- 
cause of the approximation in Eq. (21). Still we obtain a 
reasonable matching of the two approximation schemes 
for intermediate densities around ip = 0.64. It would be 
easy, in principle, to reconcile the two approximations 
at all volume fractions, including the crossover regime, 
by avoiding the low temperature approximation made in 
Eq. (21), but these computations would require a much 
heavier numerical treatment. 



B. Relaxation time and glass fragility 

The calculation of the liquid complexity, apart from 
signalling the glass transition, can be related to the re- 
laxation time of the liquid by making use of the Adam- 
Gibbs scenario [72-74] , which relates the relaxation time 
r of a liquid to its configurational entropy. In the case of 
our mean-field vision of the phase space of the liquid, the 
configurational entropy can be identified with the equi- 




FIG. 6: (Color online) A thermodynamic Angell plot of 
Ej/Eeq, with Eg = 1 as a definition of the glass transition, 
plotted against Tg/T, for different volume fractions. The ther- 
modynamic fragility is predicted to increase rapidly with vol- 
ume fraction, in agreement with numerical simulations [39]. 



librium complexity Sgq. 
relation yields: 



In that case the Adam-Gibbs 



T(r, If) = Aexp 



T^eqiT, if) 



(31) 



where c is a constant. We can push the correspondence 
between complexity and relaxation time a bit further, 
and use it to extract informations on the fragility of the 
glass. The fragility of a glass quantifies how different 
the relaxation processes are from usual Arrhenius relax- 
ation processes, and can be observed by representing the 
logarithm of the relaxation time against the inverse tem- 
perature. A linear curve is then the signature of an Ar- 
rhenius relaxation, while steeper curves indicate a greater 
fragility, i.e. a faster increase of the relaxation time. 

In Refs. [39, 40], Berthicr and Witten showed that, 
for harmonic spheres, compressing the system leads to 
a significant increase of the fragility. Using our results 
for the complexity, we are able to check qualitatively 
whether replica theory can reproduce this feature. In- 
deed, thanks to Eq. (31), the fragility can be extracted 
equivalently from the relaxation time or from the com- 
plexity. Since the fragility is usually evaluated at the 
conventionally defined laboratory glass transition, we ar- 
bitrarily define the glass transition temperature Tg{(p) as 
the temperature at which the equilibrium complexity is 
equal to one, 12{Tg{ip),ip) = 1, which is a typical value 
of the configurational entropy at the glass transition in 
most numerical simulations and experiments (its precise 
value is immaterial for our purposes) . Using these values 
of Tg and Sg = 1, we can construct an Angell plot for 
the complexity following Ref. [75]. In Fig. 6 we show the 
inverse of the complexity, linked to the logarithm of the 
relaxation time by Eq. (31), plotted against Tg/T, for 
several densities. The fragility is the slope of the curves 
in Tg/T = 1 [75]. One can clearly see that increasing 
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FIG. 7: (Color online) Optimal number of replicas m* as a 
function of T, for several volume fractions. The inset shows 
the T — >■ limit of m* as a function of the volume fraction. 



the density drastically increases the fragility of the glass- 
former. 



C. Free-energy of the glass 

We turn now to the calculation of the free entropy of 
the glass, Sgiass- We have seen in Sec. II that in order 
to compute Sgiass {T, f) , one needs to optimize the repli- 
cated free entropy S{m,A;T,ip)/m defined in Eq. (22) 
with respect to A first, via Eq. (24), then with respect 
to 771, via Eq. (13). Calling A*{T,if) and m*{T,if) the 
optimal values of A and m, we obtain: 



Sglass{T, if) 



Sim*{T,ip),A*iT,cp);T,cp) 
m*{T,ip) 



(32) 



Starting from Eq. (32), we can deduce all quantities of 
interest, for instance the pressure and energy of the glass 
and its specific heat. The derivations are reported in Ap- 
pendix B 4. In particular, our approximations imply that 
the energy of the glass Ugiass can be computed directly 
from the knowledge of the effective potential 0e/ / via the 
following relation: 



Ugiass{T,v) = l'^Vyug{v) j drr\l 



(33) 



Similarly, we find that the reduced pressure (or "com- 
pressibility factor" ) p = (iP/ p is given by 



1 



Pglass{T,ip) = —piiq{T/m*,ip) 



m* 



G{m*,A*-T). 



(34) 



In Fig. 7, we show the behavior of the optimal num- 
ber of replicas m* as a function of the temperature, for 
different volume fractions. We see that the T — > limit 
of TO* (T, Lp) converges when Lp is not too large to a finite 
value, which we call m*jjg{(f), and which is shown in the 
inset. The replica parameter of hard spheres vanishes at 
a density ^pgcp-, the glass close packing. This point is 
the equivalent, in our mean-field picture, of the jamming 
point of harmonic spheres. The behavior of the cage ra- 
dius A* is similar to that of to* . We will study in greater 
detail the behavior of the system around fccp h^ the 
next section. 

From the knowledge of m* and A* , we can deduce the 
energy and specific heat of the glass. We show in Fig. 8 
the temperature evolution of the specific heat for three 
densities, one below (facp, one very close to it, and one 
above. The temperature dependence of the specific heat 
is qualitatively in good agreement with the numerical 
data shown in Ref. [39] , the only difference being that the 
abrupt jump predicted at Tk using the theory is replaced 
by a smoother evolution and a small maximum at the 
numerical glass temperature. 

Upon crossing the ideal glass transition, the specific 
heat undergoes a finite jump. We find that the ampli- 
tude of this jump increases continuously with ip from 
ifK- A qualitatively similar result was obtained in nu- 
merical simulations [39], where the jump of specific heat 
was studied at the numerical glass transition tempera- 
ture. The behavior of the specific heat correlates well 
with the evolution of the thermodynamic fragility dis- 
cussed in Fig. 6. Finally, we note that the T limit 
of the specific heat jumps discontinuously from to 3/2 
£^t V'GCPj which reveals that the ground state proper- 
ties of the glass phase change abruptly at the jamming 
transition, as we now study in more detail. 



IV. JAMMING POINT OF HARMONIC 
SPHERES 



We now turn to the study of the region of the phase 
diagram deep in the glass phase, inside the line TK{yy) 
and close to the jamming point (pacp, see Fig. 2. In this 
region, m*{T, ip) is very small, as we discussed in the last 
section (see Fig. 7). In principle, one could just compute 
TO* (T, ip) and then take the limit T — >■ and ip — > ipcc p 
("jamming limit"). However, both numerically and an- 
alytically, it is much more convenient to exchange the 
optimization with respect to to and A with the jamming 
limit, take the latter first, and then optimize the free 
energy. This is because, in the jamming limit, many 
of the integrals that appear in the free energy simplify 
considerably. However, the scaling of to and A in the 
jamming limit is different depending on the order of the 
limits, which emphasizes the asymmetry that exists be- 
tween both sides of the jamming point. In this section 
we discuss how to exchange the jamming limit with the 
optimization procedure. 
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FIG. 8: (Color online) Specific heat Cv of the system plot- 
ted against T/Tk, for several volume fractions. The curves 
for ip = 0.620, 0.633 are obtained with the effective potential 
method developed in this paper. The curve for — 0.650 is 
obtained with the Mezard-Parisi small cage expansion [31], 
since at this density TKif) is too high and the effective po- 
tential method is not reliable. Inset: Specific heat jump at 
Tk\ the black dots are obtained with the effective poten- 
tial approximation, the gray triangles are obtained with the 
Mezard-Parisi small cage expansion. 



A. Zero temperature, below close packing 

We consider first the limit T — > for fixed < V'GCP- 
In this case, we are bound to recover tlie results of [33] 
for hard spheres. Thus, we expect the optimal number of 
replicas m* {T, ip) and cage radius A* (T, tp) to tend to fi- 
nite values 'm*jjg{tf) and A*^g{if) when T — > 0. Therefore, 
to recover hard spheres results, one has to take the limit 
T — at fixed A and to, and optimize the resulting free 
entropy over A and m. The corresponding expression for 
the replicated free entropy is given in Appendix C and 
coincides with that of [33]. We refer to this paper for 
further details on the hard sphere case. 

We solved numerically the optimization equations; the 
result for m*j^g{if) is plotted in the inset of Fig. 7. Ap- 
proaching ipGCP by the left, tends to zero linearly: 



'■HS 



11 {ip. 



GCP 



(35) 



and A*fjg also goes to zero with to as am. Thus the close 
packing limit for hard spheres can be computed by taking 
first the limit T — > 0, and then m — with A = am. 
Optimization on A will now become an optimization on 
a. From the inset of Fig. 7 we find 



Vgcp = 0.633353, 



(36) 



and Jl = 20.7487, see Appendix E for more details. 
Therefore, we find that the location of the jamming tran- 
sition is different from the density of the hard sphere 
glass transition, which confirms the distinct nature of 
both phenomena. 
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FIG. 9: (Color online) Complexity of the energy minima 
Eo(e, (f) as a function of energy for several densities. 



In this limit (T ^ at finite to and A), it is easy to see 
that 4>ef / has an hard-core of the same size of the original 
potential. Therefore, from Eq. (33), it is straightforward 
to see that the energy of the glass vanishes: the system 
becomes indeed a system of hard spheres. From Eq. (34), 
since both piiq{T = 0,(p) and G{m,A;T = 0) arc finite, 
we find that the reduced pressure goes to a finite value 
Pgiass{T = 0, which is the reduced pressure of the hard 
sphere glass. The latter diverges at ipccp since to* — > 
at that point. One can show that for tp — > Vgcp^ since 
TO* vanishes linearly, one obtains [33] 



Pglass{T = 0, ((S) 



3.03430 yccp 
Pgcp - V 



(37) 



where the prefactor is only 1% different from the correct 
value which is the space dimension d = 3, predicted by 
free volume theory [76] and by the small cage expansion 
ofRef. [33]. 



B. Zero temperature, above close packing 

For Lp > ipGGPi the harmonic spheres have a finite 
energy at T = 0, due to overlaps. Then, because of the 
repulsion the spheres are jammed, and A* vanishes at 
T = 0. It turns out from the numerical calculation of to* 
and A*, both in the small cage expansion [31] and in the 
potential approximation (see the previous section) that 
the optimal number of replicas to* tends to zero when 
T — > 0. One finds that to = T/t and A = ma with 
constants a and r. Therefore for ip > (pccp one has to 
take the T — > limit with a ~ A/m and r = T/to kept 
constant. Optimization over to and A is replaced by an 
optimization over a and r. 

The replicated free entropy (22) simplifies considerably 
in this limit. The details of the calculations are discussed 
in Appendix D. The expression for iSo(a,T;(/9) in this 
limit is given in Eq. (D4). We first calculate the optimal 
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cage radius a with Eq. (D5), and replace the solution 
a*{(p) in Eq. (D4). We obtain the replicated free entropy 
as a function of T, iSo(t; lys) ~ SQ{a* ,t; ip). To understand 
the physical meaning of this quantity, we have to go back 
to Eq. (10) and take the hmit T ^ with m = T/t. 
Since the free energy /* reduces to the energy e* for 
T = 0, we get 

5o(r; if) = lim S{T/t; T, ip) = ^o{e*) - e* /t . (38) 

T— i-O 

Following the reasoning that led to Eq. (12), wc get 

dirSo) 



So(t, (^) 
e{T,ip) = 



dr 
^dSo 
' dT ' 



(39) 



and therefore 



eGs(v) = Ugiass{T = 0, 93) = - niin[r5o(a* 

T 

= -T*{ip)So{a*,T*{p);ip), 



(40) 



since the ideal glass energy coincides with the amor- 
phous ground state energy. We obtain the scaling of r* 
and of the energy around ipccp by means of a devel- 
opment around (pccp- The calculation is presented in 
Appendix D, and the result is: 



T*{ip) - t((^- (^GCP), 

eGs(</5) ^ e(v - VGCpf, 



(41) 



for (fi — > ipQfjp. so that the jamming limit corresponds 
to T — J- 0, as expected. We get ipccp = 0.633353, 
T = 0.00835535 and e = 0.263017 (see Appendix E). In 
addition to the scalings in Eq. (41), we obtain the follow- 
ing scaling for the complexity against the ground-state 
energy for small e: 



Eo(e,</?) 



(42) 



The function I]o(e, 1^9) is reported in Fig. 9. The energy 
of the ground state ecsi^) corresponds to the vanishing 
of the complexity. We know that TiQ^^ip) ~ <pgcp — 
ip. Since the coefficient A is constant around ^gcp, one 
obtains the quadratic scaling of the energy. 

Finally, from Eq. (34) we see that the pressure 
Pnir,« JT = 0, (p) = iimx-foipTpgiass) is proportional to 



Jass{T = 0,ip) 

T* and therefore 



P{ip) P{ip - ifiGCp) 



(43) 



with P = 0.403001 (see Appendix E). This resuh is 
indeed consistent with the exact relation Pgiass{T ~ 

0,¥')=^^^. 

In this section we have thus derived the scaling laws 
(41) and (43) first observed numerically above the jam- 
ming transition at zero temperature [38, 60], which indi- 
cate that solidity emerges continuously at the jamming 
density. These results also confirm the correspondance 
between the jamming transition observed numerically 
and the glass close packing density defined within our 
theoretical approach. 
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FIG. 10: (Color online) Scaling of m'{T,ip)/VT as a func- 
tion of the rescaled temperature \ip — ^gcp\I ^JT for = 
0.58, 0.59, . . . , 0.70. The asymptotic forms corresponding to 
Eq. (45) are also reported. 



C. Scaling around jamming 

We have shown in Sees. IV A and IV B that all the 
thermodynamic quantities are singular at kpQc p and T = 
0; for instance m* is finite below 'Pgcp while it vanishes 
proportionally to T above ^gcp- The reduced pressure 
is also finite below fGCP, while it diverges at ifccp E^nd 
is formally infinite above fGCP- Both the energy and the 
pressure vanish below (fccp while they are finite above 

•i^GCP- 

From this observation, and since all quantities are an- 
alytic at finite T, it follows that all these quantities must 
satisfy scaling relations if T is small enough and ip is close 
enough to pgcp [77]. We discuss explicitly the case of 
TO* for which we assume a scaling relation of the form 



'{T, p) = T^ 



m± 



If - pgcp] 



(44) 



where the two scaling functions correspond to the two 
sides of the transition. In the hard sphere limit T — > 
and (p < ipGGP, to recover Eq. (35) we need fh^(x — > 
00) — fix and ^ = V. For p > pccp instead, to* ~ 
T/T*{ip) and to recover Eq. (41) we need TO+(a; — )■ 00) = 
1/(tx) and — v — \/2. Finally, we are led to the 
scaling 



TO* (T, p) = vT fh± 
rh^{x — > 00) = Jlx , 
m+(a: — > 00) = :z— ■ 

TX 



\V - 'Pgcp\ 



T 



(45) 



This scaling is successfully tested in Fig. 10. 

Qualitatively similar scaling forms (with different ex- 
ponents) apply to thermodynamic quantities such as the 
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energy, the pressure, or the complexity. Although we 
do not show the corresponding scaling plots explicitly, 
the energy and the pressure are plotted and compared 
with numerical data in the following (Fig. 12). Note fi- 
nally that a scaling plot for the complexity is related, in 
the spirit of the Adam-Gibbs relation in Eq. (31), to the 
scaling plot of the numerically measured relaxation time 
reported in Ref. [39, 40]. 



V. THE PAIR CORRELATION FUNCTION 

The pair correlation function is a key quantity that 
gives a lot of insight into the microstructure of dense 
fluids or packings. In this section we derive an expression 
for the correlation function of the glass and use it to 
discuss the scaling of the contact peak of g{r) close to 
the jamming transition. 



A. General expression of the correlation function 
of the glass 



It has been shown in [33, Eqs.(67)-(68)] that one can 
compute the correlation function of the glass starting 
from the replicated free energy. This is done by adding a 
small perturbation v(r) to the potential of replica 1 and 
taking the derivative of the free energy with respect to 
v{r). The glass correlation function is then calculated as: 



9giass[r) pf^Svir)- 



(46) 



We start from Eq. (19) and we use approximation (21). 
Since the perturbation acts only on replica 1, we have 



- I3rn(f){r) — I3v{r) 



+ ^2/.f/(^) /rfVe-'^™^W-^''MQ(r), (47) 



5(m, A\T, If) = Sh{m, A) + Suq 

3(y5 
TT 

and Q{r) depends only on the potentials of replicas 2 to 
TO, so it is independent of v. Taking the derivative with 
respect to v{r) we have, using again Eq. (21): 

9giass{r) = gu,{T/m, V- r) + y?,' {v)e- ^^^'^^^^ Q{t) (48) 



Vi 



HS 



((^)e-'3"^M[l + g(r)] 



= <^(^)e-^*^^/«. 

Therefore the correlation function of the glass is directly 
related to the effective potential within the low temper- 
ature approximation we used to compute the thermody- 
namics of the system. Of course, an improved expression 
for this correlation could be obtained by using the HNC 
approximation to describe the effective liquid. In that 
case, ggiass{r) would be the HNC correlation of the ef- 
fective liquid. Still, the much simpler expression (48) is 
enough to capture the scaling of ggiassii^) around jam- 
ming, at least close to the first peak. In the rest of this 



section we discuss the scaling form that is obtained start- 
ing from Eq. (48). 

A very interesting quantity related to g{r) is the num- 
ber of contacts, defined as the number overlaps per par- 
ticle. Within our low temperature approximation it is 
given by: 



z{T,ip) = 24(p / drr ggiassir) 
Jo 

= 24ipylfg^{ip) [ drr^e-'^'^^^^M . (49) 



Using the definition of the effective potential, Eq. (15), 
we can rewrite this in a simpler way: 

/>oo 

z(T, ^) = 24^y,f/(^) / du u\{A, T; u)"-ig(A, T; u) , 

(50) 

with a function q{A, T; u) which is very similar to 
q{A, T; u); its form and the details of the derivation can 
be found in Appendix B 5. 

The starting point of the analysis of the correlation 
function of the glass are Eq. (48) and Eq. (17) where the 
effective potential is defined. In the following we keep 
the notation e~^'^'ff since this function has a finite limit 
for T — > even if /3 is formally infinite. 



B. Scaling of the peak of the pair correlation at 
finite temperature 

Here we discuss the effect of thermal fluctuations on 
the maximum of ggiass{f) near the the T = jamming 
transition at (facp- This situation was studied numeri- 
cally and experimentally in Refs. [15, 41]. We focus on 
the scaling of this maximum in a region of small T and 
for ip ~ fGCP- Fig- 11 we show the behavior of the 
maximum of g glass i^), which we call gmax, as a function 
of the density, for temperatures ranging from 10^^ to 0. 
This figure demonstrates that we are able to derive an- 
alytically all behaviors reported in Refs. [41, 43, 44, 78]. 
Namely, the density at which the pair correlation func- 
tion reaches its maximum shifts towards higher values 
when the temperature departs from and increases as 
Vt, while the maximum gmax diverges as \ip — (pGCp\~^ 
on both sides of the transition. 

As discussed in the previous section, the zero temper- 
ature limit is very different for tp < ipacp or ip > ipccp- 
Thus we now discuss these two limits separately. 



C. Zero temperature below jamming: hard spheres 

In this case we send T — > before taking the jamming 
limit TO — >■ and A = am, corresponding to hard spheres 
approaching jamming from below. Taking the limit T — 
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FIG. 11: (Color online) Evolution of the maximum of the 
glass pair correlation function with T and ip. While Qmax di- 
verges on both sides of the transition at T = as Qmax ~ 
— V^gcpI"^, this divergence becomes a smooth maximum 
at finite T near the transition, whose position shifts with tem- 
perature. This behavior compares very well with numerical 
[41, 43, 44] and experimental observations [41, 78]. 



wc get the expression of [33] 
1 



r-1 



X du u 
Jo 



qiA;uy 



(51) 



with q{A; r) defined in (CI). In the jamming limit m — )■ 
and A = am one can show (see Appendix C 3) that this 
has the scaling [33] in the region r = 1+: 



1 



r-1 



(52) 



where the function Aq, is given in Eq. (C16) and is very 
well approximated by Ao(A) = 1 — ^/nXe^'{l — erf(A)). 
For if — )■ V'gcp '^^^ m ^ Jl {(pccp^^) a-nd A* = a*m 
with a finite a* = a* {ipccp)- then the peak of the cor- 
relation function for r = 1+ has a simple scaling form: 



gglass{r) 



^GCP 



(53) 

the numerical values of the parameters can be found 
in Appendix E. The height of the peak scales as 
ggiass{r = l+,<^ ^ Vgcp) ^ 1 .09922/(.^GCP - and 
its width as ipccp ~ ^- One can obtain scaling functions: 



gglassir) {(PGCP -</') = / 



with /(A) = 




(54) 
(55) 



or alternatively: 



=f[^r- 1)(1 -f 4^g,,_(l+))] , (56) 

yglassy^ ) 



with /(A) - A„. 



(57) 



The number of contacts is the integral of the peak at 
r = which is given by 



z ^ 24ip I drggiass{r) 

) peak 
,.HS/ 



^2A^yil^{^)./Aa* dAA„.(A) (58) 
Jo 

Replacing the values given in Appendix E, we get 
z = 6.13720 for (y3 = fccp- which is slightly above the 
expected z = 2d ~ 6. 

Note that it has been shown in [33] that within a sys- 
tematic expansion in ^ya, one gets z = 2d at first order. 
Here instead, in order to match with the soft sphere com- 
putation, we are using an approximation scheme which is 
not a consistent expansion in y/a. This could explain the 
fact that z ^ 2d. It would be interesting to have a full 
computation of the next term in the small a expansion to 
check whether the equality z = 2d survives at this order. 
Unfortunately this seems to be a quite hard task. 



D. Zero temperature above jamming 

For ip > ipGCP, the appropriate limit is to send T — > 
with m ~ T/t and A = am. Starting from Eq. (D13) and 
taking this limit, the effective potential takes the follow- 
ing form (the calculation is presented in Appendix D3): 



^ 0(^r - l) + e{l-r)u I r- 

2 



Aa 



T + Aa, 

(l + 4a/r)(l + l^(^ye-^('^-i)\ (59) 



Again, since t ^ Lp — pccp, close to pacp this func- 
tion develops a delta peak close to r = 1 with height 
{p — Pgcp)~^ and width ip^pccp- The main difference 
wc observe with the hard sphere limit is that, despite the 
fact that both cases give rise to a delta function at r = 1, 
the mechanism is different since for hard spheres the sup- 
port is concentrated on r = 1+ while for soft spheres it 
is concentrated on r = 1^. Another interesting obser- 
vation is that the theory predicts that ggiassif) should 
vanish exactly for r < ; this means that no pair of 
spheres can have an overlap larger than ^^^^ ■ Indeed 
in the limit r — >■ we recover the GCP at which no 
overlaps are present. Unfortunately, a direct numerical 
test of this predicion is impossible: this is because when 
1 — r = the exponential term in Eq. (59) is equal 

to exp(— 1/(t -I- 4a)). Since both r and a are already 
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extremely small (of the order of 10~^ at best), ggiassir) 
is extremely small in the region where we predict it to 
vanish, and numerically the signal to noise ratio becomes 
too large in that region. 

Again in the limit ip — > fccp obtain the scaling 
forms 



9glass{r) (if - ipGCp) = / 
,HS 



1 



/(A) = 



f - fGCP 
2 



1 - 



4a*A 



or alternatively: 




)2 . Q 



(60) 
(61) 

(62) 
(63) 



The number of contacts is easily obtained from Eq. (49) 
and Eq. (59) by a change of variable: 

zo(^) = z{T = 0,^) = 2A^y"^{^) /'drr^e-'^^^^M 



(64) 



Jo 



with Gq^ (a) given in Eq. (C7). The expression above is 
linear in r at small r. Therefore, for (p — > (^^pp. since 
T {(fi - (pGCp), we have zo{ip) = zq{lpgcp) + C{(p - 
Vgcp)- This result is at odds with the numerical finding 
of a square root behavior [43, 44]. We comment further 
on this discrepancy in Sec. VII. At ipGCPi we have 

zoifGCp) - 8^2/f/(^)G^^(a + r/4) = 6.13720 , (65) 

which is differs slightly from the expected isostatic value 
z = 2d, for the technical reasons already discussed above 
in the hard sphere case. 

Note that this value of zo{(pGCp) is identical to the 
one obtained on the hard sphere side (as it can be easily 
proven analytically) , despite the fact that the two defini- 
tions of z are not equivalent. In the soft sphere case we 
defined z as the number of overlaps. On the hard sphere 
side, this number is strictly zero, therefore we defined z 
as the integral of the delta peak at r = 1+. 



VI. SCALING AROUND JAMMING: 
COMPARISON WITH NUMERICAL DATA 

In this section, we compare the prediction of the theory 
with numerical data, in particular to test the scaling in 
temperature and in Sip ^ (p — ipGCP ai'ound the glass 
close packing point. 



A. Details of the numerical procedure 

We investigated a system of iV = 8000 identical har- 
monic spheres. The system was prepared using the stan- 
dard procedure of Ref. [38] . We started from a random 
configuration at high density; we then minimized the en- 
ergy, and then reduced slowly the density, at each step 
minimizing the energy again, until a jammed configu- 
ration of zero energy and volume fraction ipj was found. 
This jammed configuration corresponds, according to the 
mean-field interpretation of [33], to the T = 0, p — > oo 
limit of a given hard sphere glass, or equivalently to the 
T — > 0, P — )- limit of a corresponding soft sphere glass. 
Starting from that configuration, we prepared configura- 
tions at different T and p by using molecular dynamics 
simulations as in Ref. [40] to thermalize the system inside 
the glass state selected by the initial jammed configura- 
tion. Since we always used very small T and ip ipj, 
the system is not able to escape from that glass state, so 
it always remains in the vicinity of the original jammed 
configuration. This is a crucial requirement to numeri- 
cally obtain scaling behavior near ip>j as microscopic re- 
arrangments would directly affect the location of ipj in 
the simulations [79]. We tested the absence of such an 
effect in our simulations by repeating the energy mini- 
mization starting from some equilibrated configurations 
at different T and ip, and checking that we always found 
the same value of ipj (within numerical errors). For the 
particular series of run described below, we estimated 
ipj = 0.643152 ± 0.000020 by fitting the zero tempera- 
ture energy using e{p) (x {ip — ipj)'^- 



B. Difficulties in the comparison with the theory 

When comparing the theory with the numerical data, 
we face two difficulties. First of all, the value of ipj 
depends on the numerical protocol and on its particu- 
lar realization we used, and therefore cannot be directly 
compared with ipGCP- Note that (pGCP should represent 
an upper bound for ipj, whatever the protocol, but the 
value fGCP = 0.633353 we report seems to contradict 
this statement. The reason is simply that our theoretical 
computations are based on the HNC equation of state for 
liquid hard spheres, which (as discussed above), underes- 
timates the value of ipGGP'- a- better result is obtained us- 
ing the Carnahan-Starling equation of state, which gives 
'fiGCP — 0.683; consistently with the value of ipj found 
above [33]. To avoid confusion, here we stick to the use of 
the HNC equation of state; the theoretical calculations 
can be easily repeated for any other equation of state, 
and the scaling around ip>GGP is unaffected by this choice. 
Once again we stress that the absolute values of density 
reported in this paper should always be taken with this 
caveat in mind. Since we are interested here mainly in 
the relative distance to jamming, which is ip — ipj for the 
numerical data, and ip — ipGCP for the theory, the abso- 
lute values of density are not crucial for our purposes. 
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The second difficulty is the following. For hard spheres, 
the theory predicts that the reduced pressure diverges 

as Pglass ^ SfGCp/i^GCP " <^), whilc Vgi^ssi'P) 

l.l/{ipGCP — f) for 'P — > ^~GCP^ ^'^'^ Appendix E. As 
one can easily check, these scalings imply that the the- 
ory violates the exact thermodynamic relation p^f^^ = 
1 +4iy9y^£,^((y9), as already noticed in [33]. In particular, 
the divergence of p^f^^, is correct, while the coefficient of 

the divergence of y^i^gs iv) is overestimated by a factor of 
^ 1.46562. We will comment on the origin of this prob- 
lem in Sec. VII. This error affects the prefactors in the 
scalings of ggiass{f) and of z around facp- To correct 
for this error, we rescaled the 5<p of the numerical data in 
such a way that the coefficient of the divergence of y^i^^g 
is the same as the theoretical one. 

To summarize, when comparing numerical data with 
the theory, we consider Sip = ip — facp for the theory, 
and 5lp = 1.46562 (1^9 — Lpj) for the numerical data such 
that the values of y^£jj {Sip) for theory and simulations 
are exactly the same close to Sip = Q. Note that this 
adjustment is based solely on the analysis of the hard 
sphere side of the jamming transition {5lp < 0), but we 
find that it allows to describe the full scaling also for 
positive Sip. 

If the reader is uncomfortable with the above reason- 
ings, another way of thinking is that we need to ad- 
just two free parameters (the position of the transition 
and the amplitude of the rescaling of 5(p) to obtain the 
best comparison of theory and numerical simulations. 
While the location is ipj was always adjusted in previ- 
ous work [32, 38, 43, 44], we need to adjust also one 
prefactor as we seek to compare theoretical predictions 
to simulations in absolute values, and not only at the 
level of leading diverging contributions. 



C. Energy, pressure, average coordination 

We start our discussion by the simplest observables. In 
Fig. 12 we plot the average energy Ugiass{T,(p) and the 
inverse reduced pressure l/pgiass(T, ip) as functions of6ip 
for several temperatures. The agreement between theory 
and numerical data, with the rescaling of Sp discussed 
above, is nearly perfect. The scaling around jamming is 
clearly visible in the figures. For instance, UgiassiT, ip) 
tends to a finite value for Sip > 0, while for Sip < it 
goes to zero as a power law, since in this case the system 
becomes a hard sphere glass. Similarly, the reduced pres- 
sure is finite for 5p < 0, while it diverges proportionally 
to /3 for Sp > 0, since in this case the pressure is finite 
at zero temperature. At finite temperature, the curves 
interpolate between the two regimes. Scaling functions 
similarly to the ones shown in Fig. 10 for m* can easily 
be constructed. 

In Fig. 13 we report the average coordination num- 
ber z(T, as a function of Sp for several temperatures. 
At T = 0, the average coordination jumps from to 6 
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FIG. 12: (Color online) Energy UgiassiT, ip) (top panel) and 
inverse reduced pressure l/pgiass{T, p) (bottom panel) as 
functions of distance from jamming Sp, for several temper- 
atures. We define Sip — ip — pacp for the theory, and 
Sip = 1. 46562 (</3 — ipj) for the numerical data. 



(for numerical data) or to 6.13720 (for the theory). For 
Sp > 0, the average coordination grows linearly in Sp 
for the theory, while it grows as Sp^^'^ for the numerical 
data. Therefore, the theory fails to capture the correct 
evolution of this structural property at Sp > and T = 0. 

A more detailed description of what happens is ob- 
tained by looking to the data at finite T. When temper- 
atures become too large, e.g. at T = 10~^, the theory 
eventually fails because of the low-T approximation we 
made on the liquid theory in Eq. (21). For any fixed tem- 
perature T < 10~^, we observe that the theory describes 
perfectly the numerical data for Sp < and also for pos- 
itive Sp up to some crossover value Sph(T). Around 
Sph{T) the theory starts deviating from the numerical 
data. The data in Fig. 13 suggest that Sph{T) is an 
increasing function of T, i.e. that for higher tempera- 
tures the theory performs better, if the temperature is 
low enough that approximation (21) makes sense. In- 
deed, we expect that for Sph{T) — > for T — 0, since 
at T = the theory fails to capture the behavior of the 
contact number at Sp > 0. We will come back to this 
point below. 
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FIG. 13: (Color online) Average number of contacts z as 
a function of distance from jamming 5ip, for several tem- 
peratures. We define Sip = ip — pacp for the theory, and 
Sip — 1.46562(iy5 — ipj) for the numerical data. 



S(fi = -0.00075 - Tft = 




0.9992 0.9996 1 1.0004 1.0008 1.0012 



r 

FIG. 14: (Color online) Pair correlation ggiass{r) just below 
jamming predicted by theory (full lines) and measured in nu- 
merical simulations (symbols). Here Sip — —0.00075. 



D. Correlation function 

We now report data for the scaling of the contact peak 
of Qgiassir) near jamming. We first consider a fixed Sip < 
0, and change the temperature. In Fig. 14 we report 
data for 6ip = —0.00078 for several temperatures. In this 
regime the theory works very well down to T = 0; this is 
expected since we know that the theory works at T = 
for hard spheres [33]. 

Next, we consider a fixed 6ip > and plot ggiassif) 
for several temperatures. This is done in Fig. 15 for four 
different values oi dip. Here, we observe that, like for the 
average coordination (Fig. 13), the theory works very weh 
at moderately large T. However, on lowering T, at some 
point we observe that the numerical ggiassif) saturates 
to a limiting curve which has a maximum at r < 1. By 
contrast, the theoretical curves slowly evolve towards the 



T = theoretical limit, which is a half-Gaussian centered 
in r = 1, with its r > 1 part removed. 

Therefore, below a certain temperature Tfi{Sip) (the 
inverse of the function Siph{T) mentioned above), a de- 
viation between theory and numerical data is observed 
close to the maximum of ggiassif)- Clearly, the deviation 
between data and theory is a smooth crossover, so deter- 
mining ThiSip) is not easy. Here we choose the following 
procedure. The numerical data accumulate on a master 
curve as T — >■ 0. Therefore we define ThiSip) as the value 
of temperature at which the theoretical curve best fits 
the r — > numerical curve. As an example, in the upper 
left panel of Fig. 15, we see that the theoretical curve 
for T = 10^^° fits very well the numerical curves for 
both T = 10^1", 10-", while the T = 10"" theoretical 
curve is quite different. We fix — 10^^'^ for this value 
of Sip. Clearly the ambiguity on the precise numerical 
value of Til is rather large. Despite this, we are able to 
qualitatively confirm that Thi6ip) is an increasing func- 
tion of Sip, as already suggested in the previous section. 
The function Th iSip) determined from the pair correlation 
functions is reported in Fig. 16. It emerges continuously 
from zero above ipccp- A comparison with the temper- 
ature scale given by the glass temperature shows that it 
is orders of magnitude smaller than Tk. Thus it is only 
in the small regime of T < Th and ip > ipccp that our 
the effective potential approach misses some important 
physics, as we discuss below in Sec. VIE. 

Finally, in Fig. 17 we show tlia-t Qgiass 

(r), at r = 0, 

follows the scaling relations predicted by the theory, 
Eqs. (54) and (60), for \5ip\ — > 0, with different scal- 
ing functions on the two sides of the transition. While 
the scahng is perfect for Sip < [33], for Sip > there 
is a small deviation around the maximum as already dis- 
cussed, which is however very small in the scaling regime 
for Sip — > 0. 



E. Discussion 

The main discrepancy between the theory and numer- 
ical data is in the region Sip > Q and very small tempera- 
ture T < ThiSip). Here, the average coordination scaling 
z (X Sip^^^ at r = is not obtained. Furthermore, the 
shape of the peak of ggiassif) is only partially captured 
by the theory. We tentatively attribute these discrep- 
ancies to the particular nature of the vibrational modes 
at r = and Sp > 0. It has been shown that in this 
regime the low-frequency spectrum is dominated by spa- 
tially correlated vibrational modes. These modes cannot 
be described by our approach, in which vibrations are 
assumed to be Gaussian and only two-body correlations 
are taken into account. Therefore it is to be expected 
that the particular features of the jamming transition 
that are related to these correlated soft modes are not 
well reproduced by our theory. 

On the hard sphere side, the soft modes induce a 
square-root divergence of the pair correlation function 
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FIG. 15: (Color online) Pair correlations ggiass{r) above jamming predicted by theory (full lines) and measured in numerical 
simulations (symbols). The agreement is very good for T > Th{Sip), the values of Th{Sip) axe reported above the figures. 
Below Th{ifi), the numerical curves do not evolve anymore, while the theoretical curves continue to evolve towards the T — 
theoretical limit. 



Qgiassir) oc 1 / T — 1. Yet, for (5(^ — > 0~, this square root 
singularity is well separated from the contact peak [33]; 
therefore the shape of the scaling function for the contact 
peak is not affected by the soft modes and it is indeed 
correctly reproduced by the theory (see [33] and Fig. 17). 

On the other side of the transition and for small 
5lp > 0, it is commonly assumed that both the con- 
tact peak and the square root singularity are shifted by 
Under this assumption, the square root contribution 
should become ggiass{,r) ^ ^I^Jr — 1 + abip. Inserting 
this form of ggiass{r) into Eq. (49) yields the result that 
z = zo + cSip^^'^, where Zq is the contribution of the con- 
tact peak. While zq is clearly unaffected by the shift, the 
scaling of z with 6(1) directly stems from the square root 
singularity of the pair correlation function. Therefore, 
the scaling of z is dominated by the soft modes contribu- 
tion, which might explain why it is not well captured by 
our theory. 

On the other hand, inserting the same expression in 
the general formula for the internal energy, Eq. (26), one 
finds that cqs ^ eSip^ + dS(p^/^, the first term being 
the contribution of the contact peak, the second being 
the one of the square root singularity. Therefore, for the 



energy the contribution of the anomalous modes is sub- 
dominant, and for this reason we observe that the theory 
reproduces well the numerical data in this case. The 
anomalous correction Sip^^'^ is of course not reproduced 
by the theory, that gives an expansion of ecs in integer 
powers of Sip, but these are subdominant terms. 

We also suspect that the shape of ggiass{r) near the 
peak is also influenced by a mixing of the contact peak 
with the square root singularity. While both contribu- 
tions are well separated for (p < Lpj, this is not necessarily 
the case above jamming, which likely explains the small 
deviations observed in Figs. 15 and 17, but we could not 
find an empirical or scaling argument to disentangle both 
contributions. 

Finally, it is interesting to remark that these deviations 
between theory and simulations are washed out by a finite 
temperature T > Th{S(f). We suspect that Th{Sip) repre- 
sents a temperature below which the many-body direct 
interactions induced by the presence of many replicas be- 
come relevant (see the discussion in section II C). Above 
this crossover temperature, vibrations are probably much 
less correlated and our approximations become correct. 

A possible interpretation is that the low-temperature 
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FIG. 16: (Color online) The crossover temperature Th{<p), as 
determined from Fig. 15, and compared with the Kauzmann 
temperature from Fig. 5. 

correlated vibrations that are relevant below Th{S(p) are 
induced by the presence of strongly correlated soft vi- 
brational modes at T = [14, 32, 38, 56-60]. Interest- 
ingly, the fact that Tfi{Scp) — >■ for 6(p — > 0+ gives a 
hint that the correlations of these modes at Sip = are 
extremely fragile, such that an arbitrary small tempera- 
ture is able to destroy them and restore the agreement 
between theory and numerical data. This observation 
is consistent with the numerical finding that around the 
jamming point, strong anharmonicity is present [80, 81] 
and the Hessian matrix of the potential does not capture 
the correct density of states of the vibrations [81]. 

VII. CONCLUSIONS 

In this final section we summarize our results and we 
discuss some perspectives for future work. 



A. Summary of results 

Understanding theoretically the physical properties 
of dense packings of soft repulsive particles is a fully 
nonequilibrium problem, because the jamming transition 
occurs in the absence of thermal fluctuations deep inside 
a glassy phase where phase space is dominated by the ex- 
istence of a large number of metastable states. Athermal 
dynamical protocols sample packings with highly non- 
uniform weights [46, 62], while thermal protocols fall out 
of equilibrium at the glass transition [43] . 

Nevertheless, we showed that it can be successfully ad- 
dressed using equilibrium statistical mechanics tools. We 
have developed a mean-field replica theory of the jam- 
ming transition of soft repulsive spheres which satisfac- 
torily derives, from first principles, the existence and lo- 



cation of a jamming transition. This transition, that 
happens only at T = 0, is exactly the same as the 
SAT/UNSAT transition of random constraint satisfac- 
tion problems [48, 49, 52]: it is the point where the Parisi 
parameter m* that characterizes the IRSB glass phase 
goes to zero, the reduced pressure of the hard sphere 
glass diverges, or equivalcntly its pressure and energy 
become finite. Within our approach, the jamming tran- 
sition occurs deep into the glassy phase, and is thus a 
phenomenon which is physically distinct from the glass 
transition itself [33, 35]. Finally, we have shown that the 
average coordination jumps from zero to ^ 6 at the tran- 
sition, and we derived scaling functions that describe well 
the contact peak of the pair correlation function. 

We want to stress here again that the absolute values 
of density reported in this paper are affected by the fact 
that we used the HNC equation of state for the liquid [33] , 
and thus quantitative agreement with simulations cannot 
be expected. Still, we found that the scaling around jam- 
ming is almost insensitive to the particular equation of 
state that is chosen for the liquid. We also investigated 
the region of small T around the jamming transition. We 
showed that the theory reproduces quite well the scaling 
of numerical data with T and 6ip, the distance from the 
jamming transition. On the hard spheres side of the tran- 
sition {S(p < 0) the theory works well down to T = 0. On 
the soft sphere side {S(p > 0) we have identified a tem- 
perature Tfi{5tp) below which, according to our interpre- 
tation, correlated vibrations [56, 58] that are neglected 
in our theory become relevant for some observables, and 
produce a series of interesting scalings that are not cap- 
tured by our theory [14, 32, 59]. The most striking of 
these is the scaling of the number of contacts which our 
theory fails to reproduce. It is a major open problem to 
try and include in the theory a correct description of the 
divergent correlations in the vibrational modes [56, 58] 
that are relevant for 6(p > and T < Th{5Lp). 

B. Technical aspects 

There are several more technical aspects on which the 
theory could be improved, even in the region where it 
works quite well. 

Our theory falls within the general framework of the 
replica method, using the molecular liquid formulation 
of [30, 31]. However, to be able to describe the jam- 
ming transition, we had to develop a new approximation 
scheme, that is based on the effective potentials formu- 
lation of [33] , and allows to write the molecular liquid as 
a simple liquid with effective potentials involving an ar- 
bitrary number of particles. Our approximation is based 
on the following steps: 

1. As in [33], we neglect the three (and more)-body 
interactions, and just keep the two-body effective 
potential. 

2. On top of this, we observe that the effective po- 
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FIG. 17: (Color online) Scaling functions at T = above (left) and below (right) the jamming transition showing the convergence 
of the first peak near r = 1 to a delta function with asymmetric scaling functions on both sides of the transition. To ease 
visuaUzation, we show the evolution of g{r)/gmax, where gmax can be read from Fig. 11. 



tential is given by the original potential (rescaled 
by a factor m) plus a small perturbation, therefore 
we use liquid perturbation theory to write the free 
energy of the replicated liquid as in Eq. (19). 

3. Finally, we performed the low temperature approx- 
imation Eq. (21), in order to simplify the numerical 
calculations. 

The last two approximations have been done mainly for 
practical reasons. However, we found that they arc re- 
lated to the first one. Indeed, if we get rid of approxi- 
mation 3, we introduce an explicit dependence of yuq on 
temperature, hence on r in the limit T — > and S^p > 0. 
It is easy to realize that this dependence affects the scal- 
ing close to jamming. An additional term ^/t appears 
in the small r expansion (see Appendix D2), which im- 
plies that r oc dtp^ and e (x Sip^ . An inspection of the 
small r expansion reveals that if yuq does not depend 
on r, then a delicate cancellation happens to eliminate 
the T^/^ term and produce the correct scalings. The cru- 
cial observation is that the anomalous term contains the 
derivative of yuq with respect to temperature, which is re- 
lated to a three-body correlation function. Therefore we 
suspect that in a full treatment this term is cancelled by 
a term coming from the three body interaction potential. 
Therefore, getting rid of the third step of approximation 
might require taking into account three-body correlations 
as well, therefore getting rid of all steps at once. 

This is of course an extremely interesting direction for 
future research, mainly because we suspect that getting 
rid of approximation 2 should allow to eliminate the vi- 
olation of the exact thermodynamic relation p^fg^ = 

1 + ^wfiassiv) that happens in the theory and forces us 
to rescale Sip in order to compare with numerical data. 
It is also possible that taking into account interactions 
involving several particles one could capture at least par- 
tially the long range correlations at jamming. 



Finally, our theory neglects the existence of rattlers 
(particles with less than 4 contacts, which are therefore 
not mechanically blocked in jammed packings) since the 
cages are assumed to be identical for all particles [33]. 
Therefore the rattlers are approximated as jammed par- 
ticles. This should not influence the shape of g{r) close 
to the contact peak since typically rattlers are not in 
contact with other particles, so they do not contribute 
in this region, which was the main focus of this work. 
The theory can be in principle extended to include cage 
heterogeneity, and therefore describe rattlers. The com- 
putations will be more involved because a distribution of 
cage radii P{A) should be included. 



C. Some general directions for future work 

Our approach is general enough that it can be sys- 
tematically improved and generalized to various models 
(e.g. Hertzian potentials, or truncated Lcnnard- Jones 
potential). Moreover, recent work reported on a method 
to compute the shear modulus within the replica ap- 
proach [82, 83]. This could allow to add the external drive 
to the theory, and obtain a complete theoretical picture of 
the jamming transition in the three dimensional temper- 
ature, density, stress jamming phase diagram originally 
proposed by Liu and Nagel [13]. 

Finally, we believe that in the long term these studies 
could open a way towards a quantitative renormalization 
group treatment of structural glasses, following recent 
studies that formulated the renormalization group for lat- 
tice models [67, 68]. If successful, this program could lead 
to a correct treatment of long-range correlations in the 
glass phase [56, 58], and in particular at the jamming 
point, and could also allow us to address the behavior of 
pair correlations over a broader range of distances, since 
pair correlation functions also exhibit singular behavior 



21 



beyond the contact peak [43, 44], while the structure fac- 
tor [84] and isothermal compressibility [85] also display 
intriguing anomalies at large scale which remain to be 
predicted theoretically. 
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Appendix A: Complexity and the free-energy of 
basins 

The replica calculation presented in this paper is based 
on the assumption that phase space can be decomposed 
into a set of "pure states", which is the basis of Eq. (2). 
The precise definition of these states in statistical me- 
chanics is however tricky and not very practical for nu- 
merical implementation. Although less rigorous, a more 
intuitive presentation in terms of the minima of the po- 
tential energy function, sometimes called "inherent struc- 
tures" [86] or "basins", is very enlightening and allows 
an explicit connection with many numerical procedures 
encountered in the litterature to study the jamming tran- 
sition [38, 69]. For the sake of completeness, we include 
this discussion below. It is worth noting, however, that 
a partitioning of phase space in basins always depends 
on the details of the way these basins are defined. We 
expect that the particular definition we use below will be 
a good approximation of the equilibrium decomposition, 
at least for small enough temperatures. 




FIG. 18: A partitioning of phase space into basins, labeled by 
a representative configuration Ro, which we assume to have 
minimal energy in the basin. Of course the shape of the basins 
is much more complex due to the high-dimensional nature of 
phase space. The above is just a very schematic representa- 
tion. 



in the basin at the end of the procedure. Note that the 
partition defined in this way is, a priori, dependent on 
the precise procedure used to perform the energy mini- 
mization or the compression. It is reasonable to expect 
that the ambiguity due to this is small and negligible in 
most cases. 

The potential energy of the harmonic spheres is V[R] = 
Si<i ^^'^ define the free energy of a basin 

as 

g-/3w/[fl„] ^ JdR GiR, Ro) e-'5^[««l (Al) 

where G{R, Rq) ~ 1 if and only if configuration R belongs 
to basin Rq. Note that for /? — > oo we have 



1. Definition of basins and their internal free 
energy 

We suppose that the phase space of all configurations 
of the spheres R can be unambiguously partitioned in 
basins. Each basin is labeled by a representative config- 
uration, Rq, see Fig. 18. For instance, one might think 
that a basin represents the set of configurations that end 
up in Rq after an energy minimization (Rq is an energy 
minimum in this case). Or, for hard spheres, we can as- 
sume that a basin is the set of configurations that end 
up in Rq after a fast compression. One can, alterna- 
tively, minimize the energy first, then if the energy is 
zero (the configuration is not jammed), perform a fast 
compression, in such a way that Rq is always a jammed 
configuration. In any case, Rq will have minimal energy 



g-/3Af/[fio] _^ 



VbiRo) = 
„-l3Ne[Ro] 



if e[Ro] = 
if e[i?o] > 



(A2) 



where Vb is the volume of the configurations of the basin 
that are compatible with the hard sphere constraint (its 
logarithm being the internal entropy of the basin). In- 
deed, if the basin contains zero-energy configurations 
(e[i?o] = 0), then at zero temperature its free energy 
is dominated by the number of those configurations. In 
the opposite case, if there are no hard sphere configura- 
tions (e[_Ro] > 0), it is very reasonable to assume that 
the minimum of the energy is unique and if we assume 
that Rq is the configuration of minimal energy inside the 
basin, then the weight is just dominated by its energy 
e[i?o] = V[ll>]/N. 
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2. The replica method and the complexity 

The replica method allows to compute the following 
object (for a given temperature T and volume fraction 
ip), see Eq. (9): 



(A3) 



-Ro 

We can then compute a configurational entropy (or com- 
plexity) E(m;T, (^) according to Eq. (12): 

d 



I](m; T, if) ~ —m^ 



dm 



1 



logZ„ 



^l3mNf[Ro] 

-H—z log 

Ro 

Ro 



-l3mNf[Ro\ 



where 



-^mAr/[_Ro] 



Pro 



(A4) 



(A5) 



Therefore, the complexity is the entropy of basins, 
weighted by the weight Pj^^^- The weight is constructed 

by assigning an "energy" Nf[Ro] to each basin, and then 
constructing a Boltzmann weight at the effective temper- 
ature T = T/m. 

3. Special limits 

Equation (A4) allows to give simple interpretation to 
the different limits we considered in Sec. IV, as follows: 

1. Consider m = 1. Then the basins are weighted 
by their equilibrium free energy at temperature T. 
The associated complexity 'E,eq(T,ip) = S(m = 
l;T,ip) corresponds to the number of basins that 
dominate the equilibrium partition sum. In other 
words, if we take an equilibrium configuration at 
given (T, ip), this configuration belongs, with prob- 
ability 1 for oo, to one among exp[NJ^eq] 
basins (all the other basins having a negligible prob- 
ability at equilibrium). This equilibrium complex- 
ity is reported in Fig. 4. 

2. Consider taking /? — > oo first, then m = 1. In 
this case we get, at low enough density where hard 
sphere configurations exist (omitting the normal- 
ization) 

Pjig oc Vb[Ra] . (A6) 

The basins are weighted by their hard sphere vol- 
ume. Therefore Sf^^(^) = I](m = 1;T = 0,ip) 
is the equilibrium complexity of hard spheres: if 
one takes one hard sphere configuration at equi- 
librium at volume fraction if, the latter will belong 



with probability 1 to one among exp[A^I]^'^] basins. 
This quantity is reported in Fig. 16 of Rcf. [33] and 
in the inset of Fig. 4. 

3. Consider taking m ~ 0. In this case, the basins 
are not weighted. Therefore, I](m = 0; T, ip) is just 
the logarithm of the total number of basins, of any 
energy and entropy. 

4. Consider taking /3 — > oo first, again in the region 
of ip where hard sphere configurations exist. Then, 
the weight becomes 



Pro °^ M^oY 



(A7) 



which means that basins are weighted by their 
volume of hard spheres configurations raised to 
the power m. In particular, those basins that do 
not contain any hard sphere configuration (those 
of positive energy) have zero weight and are not 
counted. Now, take m — > 0. The resulting weight 
is 



oc 0{vb[Ra] > 0} , 



(A8) 



in other words we give uniform weight to all basins 
that contain at least one hard sphere configuration. 
The corresponding complexity is 



I]^^(^)- hm lim E(m;T,^) 



(A9) 



and it counts the total number of basins for hard 
spheres. 

5. Finally, consider taking the limit /3 — >■ oo, with 
the effective temperature r = T/m held constant 
(meaning that m — > at the same time). In this 
case, the weight becomes 



Pro 



-NelRo]/T 



(AlO) 



Therefore, I]o(t, (^) = \hnT^o'S{m = T/T;T,ip) 
counts the basins weighted by their minimal en- 
ergy at temperature t, irrespective of their entropy 
(in particular basins that contain hard sphere con- 
figurations have weight 1 in this case). Obvi- 
ously, lim^^oSo(r,(p) = T.^^{ip) for ip < ifacp- 
Moreover, one can perform a Legendre transfor- 
mation with respect to r to eliminate r in favor 
of its conjugate variable (the basins minimal en- 
ergy) and therefore obtain Eo(e, (/?), which counts 
the basins of minimal energy e at density ip. Again, 
lime^oSo(e, <^) = T;^^{ip) for (p < (pccp, while 
for if > ipGCP, one finds that I]o(e, (ys) vanishes 
at a positive energy ecsif) which is the global 
minimum energy for that density and increases as 
eGs(v) oc (ip-ipGCp)'^ [54]. The quantity T,oie,ip) 
is reported in Fig. 9. A similar object was recently 
measured numerically [69]. 
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Appendix B: Some calculations and simplications 
1. The function q 



The function q{A^ T; r) that appears in the effective potential Eq. (17) is given, using the definition and introducing 
bipolar coordinates, by: 



1 



rV An A Jo 

The Gaussian integral can be explicitly computed and the result is 



du 



6 — C 



-/30(u) 



qiA,T;r) 







> - 1' 




hErf 






2s/A_ 



Erf 



r + 1 



_(r-i)^ (r + AAB) 



2r(l + 4A;3)3/2 



2y/A 
Erf 



(l + 4A/3)0Fr 
r + 4Af3 



2r(l + 4A;3)3/2 



Erf 



2^v4(l + AAP) 

r-AAp 
2^A{l+AA^ 



Erf 



Erf 



1 - r 
2^A{1 + 4A/?) 

1 + r 
2^/^(1+4/1/3) 



(Bl) 



(B2) 



2. Link with the Mezard-Parisi small cage expansion 

Here we show that our effective potential approximation reproduces the result of [31] in a small A expansion at 
fixed /?. Starting from Eq. (Bl), we note that for small A, r' is small, and we expand the potential 4){r ~ f") (which 
is assumed to be differentiable) for small r': 



- e-'^'^Wil - pAA(j){r) + p^A{y(j){r)f + ©(A^)]. 
Plugging this in (23) and using an integration by parts, one easily obtains: 



(B3) 



G(m, A; T) = 3/3A(l - m) / drr"^ e' ^""^^''^ /^(t>{r) 

Jo 



(B4) 



This result, together with Eq. (22), reproduces the result of [31], with the additional approximation we made for the 
correlation function of the liquid, Eq. (21). Note that without the latter approximation, wc would not reproduce 
the result of [31] exactly; this would require an expansion around the center of mass of the molecule and not around 
replica 1 as we did in Eq. (15). 



3. Simplification of the replicated free entropy in the low-temperature approximation 



We show here the details of the simplification of the replicated free entropy in the low-temperature approximation, 
Eq. (15). Starting from Eq. (19) and using Eq. (15), wc get 



J d\gu,{T/m,v;r)Q{r) ^ y,f/(<^) J dre-P^*<^^)Q{r) ^ yf^S (^cp)^G{m, A;T) , 



where: 



G{m,A;T)^^ J d 



(B5) 



(B6) 
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Then we observe that, using Eqs. (14) and (15), we get 

If " 

- J cPxi ■ ■ ■ dPx^cfyi ■ ■ ■ d^y,np{xi ■ ■ ■ Xm)p{vi • • • y™) H e'^'^^^"-^"). 



therefore we obtain 



/ d^Xd^Yq{A,T-X -Y)"' ^ / dVg(A,T;r)" , 



G(m, A-T)^i dr {q{A, T; r)" - e-''""^^'')] , 

^0 



which leads to the desired result, Eqs. (22) and (23). 



(B7) 



(B8) 



4. Energy and pressure of the glass 

Here we show how to compute the energy and pressure of the glass, which can be computed using the general thermo- 
dynamic relations \J = —dS/d[3 and p = PP/ p = —ipdS /dip, where p is the so-called "reduced pressure" or "compress- 
ibility factor". From expression (13), together with (22) and (23), we have that Sgiass{T, tp) = S{m* ,A*;T,(p)/m*, 
where m* and A* are determined by optimization of S{m,A]T,ip)/m. Therefore, we do not need to take explicit 
derivatives with respect to m and A. 

We get for the energy: 



Uglass{T, ip) 



1 dS{m\A*-T,Lp) 
m* dp 
dSuqiT / m* , ip) 
d{pin*) 



q[A* , T- + 0(r)e-^™*^« 



dp 



UugiT/m\ip)-l2py{f^'iip) / drr20(r)e-'^™>M - 12^y-^(^) / drr'q{A\T;ry 



_,dq{A*,T;r) 



dp 



.HSf,^\ I ^i„„2„/ A* rr. \m'-i9qiA* ,T;r) 



^-12ipy"''{p) drr'q{A\T-rY 



10 dp ' 

where we made use of Eq. (27) for the liquid energy. Now we make use of the definition of q, Eq. (Bl), to get 

dqiA,T;r) _ 1 



du u(j){u) 



6 — 6 



dp ^^/47^A Jo 

We plug this in the equation above, and we get (exchanging the names of u and r in the integral): 

1 



(B9) 



UgiassiT, ip) = 12¥.y,f/(^) / drr^^ir) e'^^^'') 



poo 

= 12(^2/;f/(^) / drr20(r)e- 
Jo 



rV AttA Jo 

04>cff{r) 



du 



g 4A* — g 4A* 



q{A*,T;u) 



m* — 1 



(BIO) 



where we used the definition (17). This is the desired result, Eq. (33). 
For the pressure, we get the simple result 



PglassiT, (p) 



ip dS{m*,A*;T,ip) 
m* dip 

—PliqiT m ,ip) 

m* m* 



(Bll) 



G{m\A\T). 
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5. Simplified expression of the number of contacts and of the energy 

Here we derive a simplified expression for tlie number of contacts. Using Eq. (15) we obtain 



Jo uViTiA Jo 

24^2/;f/(^) r du u\{A, T; uY'-^ q{A, T; u), 



dr r 



6 — C 



-I3{l-rr 



(B12) 



where we defined the function 

qiA,T;u) 



dr r 



C — 6 



(1-rf 



(B13) 



Note that q(A,T;u) is equal to the part of q(A,T;u) that vanishes in the zero temperature limit. This allows to 
compute it starting from Eq. (B2): 

qiA, T- r) = q{A, T; r) - lim q{A, T; r) 

T— 5-0 

e — e 



TT (1 + 4A^)r 



_(^_i)2^ (r + 4A/3) 



2r(l +4A^)3/2 



g- l + 4Af) ^ ^' 



2r(l + 4A/3)3/2 



Erf 



Erf 



r + 4A;9 



2^A{1 + 4.AP) 

r~^Ap 
2^A{l+AA^ 



Erf 



Erf 



1 - r 



2^A(1 + AAp) 

1 + r 
2^A(1 +4yl/fy 



(B14) 



Similarly, for the energy of the glass we get 



uy/AirA Jo 
™-i / dq{A,T;u 



dr r(l — r)^ 



e — e 



g-/3(l-r)^ (B15) 
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Appendix C: The hard sphere limit 
1. Free energy of hard spheres 



Here we discuss how Eq. (23) is simplified in the limit T — > 0, with fixed m and A, that is the appropriate one below 
fGCP- Taking the limit T — >■ at fixed A, Eq. (B2) gives back the results of Parisi and Zamponi [33, Eqs. (C4)]: 



q{A; r) = lim q{A, T; r) 2 + Erf 

T-J-O 2 



r - 1 



2^/A 



Erf 



r + 1 



2^/A 



I A\ ( (r-l)-' (r+l)"^ 
g 4"^ — g 435 



(CI) 



In this limit wc obtain from Eq. (23) the following result: 

5(m, A- ^) = 5„(m, A) + 5,f/(¥.) + 4^y,f/(y.)G(™, A) 

G(m, A) ^3 drr^ [q{A; r)"" - 0{r - 1)], 
Jo 



(C2) 



These equations coincide with the ones reported in [33, Eq. (72)]. 
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2. Jamming limit from below 



The jamming limit from the hard sphere side Lp — > fQ^p is obtained by taking to — > and A ~ am [33]. In this 
hmit, the function q{A; r) goes to 1 for ?■ > 1 and to for r < 1. Using the asymptotic expansions of the error function 
we get: 



q{A; r) 





'r-1' 




^1 + Erf 


) 




2s/A_ 





I A 1 

-e 
TT r 



A (r-lf 

(r - 1) + \ — e 1^ 

V TT 



1 1 



r 1 — r 



II -rl 



(C3) 



where the function R(x) goes to 1 for x ^ oo and is proportional to x for a; — > 0, in such a way that there is no 
singularity at r = 1. 

Based on this result, we can evaluate the integral that defines G{m, A) in Eq. (C2). For r > 1, we have q{A; r) ^ 1 
and therefore 



drr^ [q('ma; r)™ — 1] = / drr^ {'m[q{ma; ?■) — 1] + • • • } = 5*1 to^ + S21tl^^'^ + ■ ■ ■ 



(C4) 



For r < 1, we have instead 

nl 



drr q{ma; rY' 







am 

TT 



2 <''-l) 

drr e ■1° 



1 1 



r 1 — r 



-i? 



1 -r 



(C5) 



The last term in square brackets in the integral gives rise to a regular expansion in powers of m for small to. Collecting 
all these results together we get 



G(to, am) - to™/2gHS(-q,)(^ ^ _^ ^^^2 + ...) + s-^^^z _^ 5-2TO,5/2 + • • ■ 



(C6) 



Then, in the limit to ^ 0, Eqs. (C2) reduce to: 



5o«^(a; ^) = - 2 (ln(27ra) + 1) + + 4^5?///,^ ((p)G^^ (a) , 

G^'^(a) = 3^ r^e-^^dr = 3 + 2a)Erf^^^^ +2ae-a^ -4a] 



Optimization over a yields an equation for a*^g{ip): 

3 



^fzo, , dGff'^(a) 
with Jq (") = 



da 



Next we calculate the leading order in the small to expansion around the previous result. We have 
S{m,a; ^) = S^^{a; ^) + ^ [3 + iipyf^^^ {ip)G^^ (a)] to log to + 0(to, (mlogTO)^, • • • ). 



(C7) 



(C8) 



(C9) 



Therefore for the complexity we get at the leading order (recall that we do not need to add the explicit derivative 
with respect to a, since a is defined by the condition that it vanishes): 



E(TO,a;v^) = -m^ — [S{m,a;ip)/m]^S^''{a;^) - ^ [3 + 4^y;f/(^)G^^(a)] . 



TO 

2 



(CIO) 



The expression above has to be computed in a = a*jfg{ip). Therefore, Sq^ {a*{ip)] ip) = T,Q^{(p) is the complexity of 
hard-spheres at to = [33], which vanishes linearly at ip — (pacp- Since m*^g{ip) is defined by T,{m,a*ffg{ip);ip) — 0, 
we get the leading term of m*jjg{(p) for m or (p ^ ^gcp' 



3 + 'iipy"Hip)Gr{a*^siv>)y 



(Cll) 
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3. Scaling of the effective potential on approaching jamming from below 



Here we discuss the emergence of the deha peak at r = 1+ in the effective potential of hard spheres for ip 
'^GCP l"^"^]' "where A = ma and m — > as discussed above. For small A, we have 



1 



r\/ in A Jo 



du 



-q{A-u) 



rn— 1 



(C12) 



with q{A] r) given in (C3). First we note that if r — 1 ^ V^, then the integral is peaked at u > 1, the function q ^ 1 
and the integral goes to 1 too. Therefore e^^'^'f^'^^^ — > 1 for r — 1 ^ vG4. 

We have then to study the region r — 1 ~ V^- In this region one can easily show that the part of the integral at 
M > 1 gives a finite contribution. If we want to study the divergent part, we can truncate the integral at u = 1. We 
then introduce the scaling variable A = (r — l)/\/4m]4 = {r — l)/(m\/4a) > and in the integral we change variable 
to u = 1 — xy^ AA/m = 1 — x\/Aa. Using the first equation in (C3), and defining 



e(s) = -[l + Erf(s)] , 



(C13) 



we get for the divergent part 



dx (1 — XV 4a) e 



e - 



71" 1 - xy/'ia 

Making use of the asymptotic large s expansion of the error function, we have 0(— s) ^ e^* /(2y^s). and we get 



(C14) 



-/30f}"/(A) 



'7rm ,/Q 

f.l/V4a 



dx (1 - xy/^) 6-2^^-=^ 
dx (1 - xVI^) 







m 1 

TT 2a; ' V 1 - xVia. 

-1 



2.T 1 - x\/4a. 



(C15) 



-A„(A) 



where we introduced the scaling function that describes the delta peak: 



A,(A) ^ 2 



l/Vio 



dxx{l — xy 4a) e 



2 -2\x-x^ 



(C16) 



Note that at (pccp wc have a ^ 10 ^; therefore an excellent approximation of the scaling function is obtained by 
setting a = 0: 



/■OO 

Ao(A) = 2 / dxx e"2Ax-x= ^ ^ _ ^Xe^' (1 - Erf(A)) 
Jo 

that allows to recover the results of [33, Eq. (C40)]. 



(C17) 



Appendix D: The zero temperature soft sphere limit 
1. Free energy of soft spheres at T = 

When if > fGCP, the appropriate T — > limit is taken with m — T/t and A — ma for constants a and r. We 
start by looking at the asymptotics of q{am, rm; r) when m goes to 0. From the definition (B2) we find that: 



q{am, rm; r) 



r(l+4a/T)3/2 



r < 1 
r > 1 



(Dl) 
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Therefore the part of the integrals in G(m, A; T) for r > 1 vanishes in this hmit. Plugging Eq. (Dl) in Eq. (B8) we 
get: 

G{m,am;Tm) - s/ [e^T^i^ - e'-^^'^^H dr , (D2) 

m^O Jo L J 

which can be written, using the definition (C7), as: 

G{m,am;Tm) - Go(a, r) = G^^(a + t/4) - G^^(t/4) . (D3) 

m— >-0 

The free entropy is given by Eq. (22). Replacing T — rm and A ^ am and taking the m — )■ Umit, we obtain: 

So{a, r; ^) = (ln(27ra) + 1) + Su,{t, + 4^y;f/(^)Go(a, r). (D4) 

As before, optimization over a yields an equation for a* (r, ip) for each r and Lp: 

3 



Jo(a*(T, (p),t) 



Jo(a, r) = a^^(a,T). (D5) 
aa 

In this limit, the free energy is much simplified and the optimization over a and r is easily carried out. 

2. Small r expansion: jamming from above 

A very important observation is that, taking the r — >■ limit, we recover the hard-sphere value of S^^ {a; (p) given 
in Eq. (C7). Indeed we expect that r tends to zero when approaching (pccp from above. Therefore, to expand around 
the jamming point, we must perform a small r expansion of the free entropy (D4). 

We start from (D4) and make use of the expressions (28) and (C7). We note that some remarkable cancellations 
happen, that eliminate a dangerous term in y^. Keeping only the lowest orders in the expression (D4) of the free 
entropy, we get: 



So{a,T;ip) =SQ^{a;ip) +TSi{a,ip) + o(t). 



J„^^(a)-6 



a 



(D6) 



The optimal value of a is defined by : 



-j^ {a*{T, p); p) + r—^ {a*{r, p), p) = , (D7) 
da oa 

and it has the form 

a* (r, p) = a*Hs (^) + rb{p) + o(r). (D8) 

Recalling that a*fjg maximizes S^^{a;(p), therefore 8^5^"^ {a; ip) = in a = a'^g, we see that the r-dependence of 
a* does not contribute at lowest order. Thus: 

So {a*{T, p), r; p) ^ S^^ {a*Hs{^); + rS^ {a*Hs{v), ^) + o{t). (D9) 

Remember that S^^ (a* (</?); = Y,Q^{p}) is the complexity of hard-spheres at m = (see Appendix C). We define 
Si{p>) = Si (a|fg(<p), p>) and we obtain 

So (r ; ^) = (p) + tS^ (p) + o(t) . (DIO) 

Using Eq. (39), we obtain 

I]o(T,^) = I]rM+2r5i(^), ^^^^^ 
e{T,p) = r'^Siip). 

Finally, reconstructing this parametric equation defining Eq, we obtain : 



Eo(e,^) =E^^(^) + 2v/e5i(v?). (D12) 
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3. Effective potential in the zero temperature limit for soft spheres 

We start from the definition (f 7) of the effective potential and we take the limit T — > with to = T/t and A ~ am. 
A direct inspection of the different terms shows that, since A is very small, one can discard a number of terms to 
obtain: 

e-/30e//(r-) ^ , / ^^^^ j.. 1 ^ ^^g^ 

r\/ AttA Jo 

. ^fu-l\ g(^-i)^ u + AAB ^( l-u \ , , 

V\/4l/ u(l+4A/3)3/2 \^ ^4A(1 + 4yl^) J 

with Q{s) defined in Eq. (Cf3). We should now make a separate analysis for r > 1 and r < 1. In both cases, 
the integral is dominated by u '-^ r. Therefore, for r > 1 also u > I and in this case q{A,T;u) ^ f, so we get 

The more interesting case is for r < 1. We have 



"V47rA 



m — 1 

^ , u-l\ u + AAB ^ I l-u 

9 + e~ ^-^6 — 

^/ u(l + 4A/3)3/2 I ^4A(1 + 4A/3) 



, (D15) 



For r < 1, the integral over u is dominated by the decay of the function away from u = \ (for u < 1). We recall 
that 0(s — > — oo) ~ , while 0(s — >■ cxd) — >■ 1. Therefore, inside the square brackets wc can neglect the first term 
^-^^^ whose decay is faster. We get, neglecting all the terms that are small for to — >■ 



^-I^^.SS^r) ^^-p(l-r) dwue'n *° l^-^TT^ ^'^ . (D16) 

r V 47raTO Jo " + 4a/ r 

The term in the exponential proportional to Xjm dominates. Rearranging terms wc obtain: 

[Tl + 4o/T-r(l+4Q/T)]^ 

e P0e//('-)^e ™. e - / <^7;, 7/, e x+4c _1 (D17) 

r Jo y^ATwmJY+Aa/r) u + Aa/r 

= -/ duu5[u-l-{r-l){l + Aa/T)]e--^ ^ ' ' (D18) 

r Jo " + 4a/T 

4a \ (l+4a/T)[l + (r-l)(l + 4a/T)]2 _x±4c.f^_i, 



T + Aa 



e-^^"^-'' , (D19) 



since the delta function forces u — l + (r— l)(l+4Q;/r), but the integral vanishes unless this value of u is positive. 
Adding the contribution for r > 1, we get the final result: 



r -l) + e{l-r)e{r ^ ' " \ '-^ L^e^^^'-'^ (D20) 

\ T + 4a J 

Note that the above expression is discontinuous in r = 1, while at any finite T the effective potential is continuous. 
We can derive the scaling function that describes the emergence of the discontinuity for T — > as follows. We start 
from Eq. (D15) and we perform a change of variable u = \ + a;-\/4A in the integral, and define A = (r — X^jyf^A. Then 
we get in the zero temperature limit 



1 / X \ ^^^^^ 

^1 + Aa/T \ + Aa/T ) 

which provides the scaling function in the region where A is of order one, and ?' — 1 ^ \/~A. It is quite easy to check 
that the function above tends to 1 for A — > oo, while it tends to 1 + 4q!/t for A — > — oo, therefore interpolating between 
the two regimes r > 1 and r < 1. 
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Appendix E: Numerical calculations 

In this Appendix we collect all the numerical results of the various calculations presented in the paper. 

1. Hard spheres 

We start from the determination of fccP: which is done in the most precise way by looking at the point where 
S^'5((p) = 0, see Appendix C and D2. We get 

Vgcp = 0.633353 , 

. (El) 
a?^(</'GCp)= 9.72187- 10-^ 

with Vul^ifGCp) = 23.6238, and from Eq. (Cll), using G^'^ {a*us{fGCp)) = 0.0170908, we get 

m*Hs{^) = KfGCP ~^)^ 20.7487((^GCP - (E2) 
From this and from Eq. (34) we can get the scaling of the reduced pressure for ip — > y^QQp- 



1 

Pglas. 



HVGCP - 

3.03430 (/3GCP 

V'GCP - V 



PhqifGCp) - "i^GCP [ VuqiVGCp) + fGCP—^ifGCp) ] G^^ {a*Hs{fGCp)) 



HS 



(E3) 



The prcfactor 3.03430 is extremely close to the space dimension d = 3, as predicted by free volume theory; the value 
d is obtained in the small A expansion [33] . 

2. Soft spheres at zero temperature 

Close to (pGCP we have T,Q^{(p) — 62.9577((/3gcp ~ 'p) and Si{(pGCp) = 3767.51, hence we immediately deduce 
from Eqs. (Dll) that: 



T*{^) = -E^^(^)/(25i(^)) = 0.00835535(^ - ^Gcp) , 
ecsiip) = (S?^(^))V(45i((p)) = 0.263017((^ - ^acpf 

The pressure is obtained from Eq. (34): 



Pglass (T = 0, (yS) = lim [pTpgiass {T, <^)] = — 
T— yO TT 



p/r,'^(^)-4^ Uf/(^) + ^^(^) j Go(a*(^),r*(^)) 



T*(¥'), 
(E5) 



and expanding around ipccp we get 



u irr (\ \ ^^GCP 
Pglass{T = Q^ip) = 



Pliq{VGCp)~'^VGCP {Vuq {VGCP) + '^GCP—^{VGCp) \ G"^ {o'Hsi'^GCp)) 



TT 

0.403001((^ - (/3GCp) 



Clearly the same scaling could have been obtained from the exact relation Pgiass{T — 0,(p) — 5^^£gs 

3. Correlation functions 

From these values we can deduce the scaling of the peak of ggiass{f) for ip — > (p^f^pi 



(E6) 

dip 



Aa*{ip)\ _ 1.09949 ^^^^ 



r*((y9) / p- ipGCp' 
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and for ip — )■ ip^^jp: 

- ^ HS, ^ 1 , /.N 1.09922 

Note that even though the prcfactors arc almost identical numerically, we were not able to prove that they actually 
coincide. I 
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